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ABSTRACT 

The goal of this paper is to investigate in N-body/SPH hydrodynamical cluster simulations the impact of artificial viscosity on the ICM 
thermal and velocity field statistical properties. In order to properly reduce the effects of artificial viscosity a time-dependent artificial 
viscosity scheme is implemented in an SPH code in which each particle has its own viscosity parameter, whose time evolution 
is governed by the local shock conditions. The new SPH code is verified in a number of test problems with known analytical or 
numerical reference solutions and is then used to construct a large set of N-body/SPH hydrodynamical cluster simulations. These 
simulations are aimed at studying in SPH simulations the impact of artificial viscosity on the thermodynamics of the ICM and its 
velocity field statistical properties by comparing results extracted at the present epoch from runs with different artificial viscosity 
parameters, cluster dynamical states, numerical resolution and physical modeling of the gas. 

Spectral properties of the gas velocity field are investigated by measuring for the simulated clusters the velocity power spectrum 
E(k). The longitudinal component E c (k) exhibits over a limited range a Kolgomorov-like scaling oc Ar 5/3 , whilst the solenoidal power 
spectrum component E s (k) is strongly influenced by numerical resolution effects. The dependence of the spectra E{k) on dissipative 
effects is found to be significant at length scales i, 100 - 300 kpc, with viscous damping of the velocities being less pronounced in 
those runs with the lowest artificial viscosity. The turbulent energy density radial profile E nir t,(r) is strongly affected by the numerical 
viscosity scheme adopted in the simulations, with the turbulent-to-total energy density ratios being higher in the runs with the lowest 
artificial viscosity settings and lying in the range between a few percent and ~ 10%. These values are in accord with the corresponding 
ratios extracted from previous cluster simulations realized using mesh-based codes. 

The radial entropy profiles show a weak dependency on the artificial viscosity parameters of the simulations, with a small amount of 
entropy mixing present in cluster cores. At large cluster radii, the mass correction terms to the hydrostatic equilibrium equation are 
little affected by the numerical viscosity of the simulations, showing that the X-ray mass bias is already estimated well in standard 
SPH simulations. 

The results presented here indicate that in individual SPH cluster simulations at least N ~ 256 3 gas particles are necessary for a 
correct description of turbulent spectral properties over a decade in wavenumbers, whilst radial profiles of thermodynamic variables 
can be reliably obtained using N t 64 3 particles. Finally, simulations in which the gas can cool radiatively are characterized by the 
presence in the cluster inner regions of high levels of turbulence, generated by the interaction of the compact cool gas core with the 
ambient medium. These findings strongly support the viability of a turbulent heating model in which radiative losses in the core are 
compensated by heat diffusion and viscous dissipation of turbulent motion. 

Key words. Galaxies: clusters: general - X-rays: galaxies: clusters - Methods: numerical 



1. Introduction 

According to the standard hierarchical scenario the formation of 
larger structures is driven by gravity and proceeds hierarchically 
through merging and accretion of smaller size halos. Within this 
scenario galaxy clusters are the most recent and the largest viri- 
alized objects known in the universe. Their formation and evo- 
lution rate is then a strong function of the background cosmol- 
ogy, thus making galaxy clusters powerful tools for constrain- 
ing cosmological models (see Voit 2005, and references therein). 
During the gravitational collapse the gaseous component of the 
cluster is heated to virial temperatures by processes of adiabatic 
compression and shock-heating. Accordingly, at virial equilib- 
rium most of the baryons in the cluster will reside in the form of 
a hot X-ray emitting gas, which is commonly referred to as the 
intracluster medium (ICM). 

A basic feature of cluster formation occurring in this sce- 
nario is the large bulk flow motions (~ 1 000 km s -1 ) induced 
in the ICM by major merging and gas accretion. The relative 



motion between the flow and the ambient gas generates, at the 
interface, hydrodynamical instabilities which lead to the devel- 
opment of large eddies and to the injection of turbulence in the 
ICM. These eddies will in turn form smaller eddies, thereby 
transferring some of the merger energy to smaller scales and 
generating a turbulent velocity field with a spectrum expected to 
be close to a Kolgomorov spectrum. Additional processes which 
can stir the ICM are galactic motions and AGN outflows, al- 
though the contribution to random gas motions from the latter is 
expected to be relevant only in inner regions of the cluster. 

Theoretical estimates for the amount of turbulence present in 
the ICM are dominated by uncertainties in the determination of 
the kinematic viscosity v of the medium. In the absence of mag- 
netic fields, classical values lie in the range Re = UL/v ~ 100 
(Sunyaev, Norman & Bryan 2003), where U is the characteristic 
injection scale and V is the characteristic velocity. In a magne- 
tized plasma Reynolds numbers are expected to be much higher 
because of the reduction in the transport coefficients and the sub- 
sequent suppression of viscosity due to the presence of magnetic 
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fields (Iapichino & Niemeyer 2008). These uncertainties in the 
estimate of the Reynolds numbers of the ICM indicate that the 
classical value can be considered as a lower limit to the level 
of turbulence present in the ICM. A conservative assumption is 
therefore to consider the ICM as moderately turbulent. 

On the observational side, turbulence in the ICM could be 
directly detected using high resolution X-ray spectroscopy to 
measure turbulent velocities through emission line broadening. 
Unfortunately this approach is still below the current limit of 
detectability. Nevertheless, indirect evidence for the presence of 
turbulence in the ICM has been provided by a number of authors. 
From spatially resolved gas pressure maps of the Coma cluster 
Schuecker et al. (2004) measured a fluctuation spectrum consis- 
tent with the presence of turbulence. Other observations which 
suggest the presence of turbulent motion are the lack of resonant 
scattering from the He-like iron K a line at 6.7 kev (Churazov et 
al. 2004) and the spreading of metals through the ICM (Rebusco 
etal. 2006). 

ICM properties are expected to be significantly affected by 
turbulence in a variety of ways. Primarily, the energy of the 
mergers will be redistributed through the cluster volume by the 
decay of large scale eddies with a turnover time of the order of 
a few Gigayears. Turbulence generated by substructure motion 
in the ICM has been proposed as a heating source to solve the 
'cooling flow' present in cluster cores (Fujita, Suzuki & Wada 
2004), the heating mechanism being the shock dissipation of the 
acoustic waves generated by turbulence. Moreover, the useful- 
ness of clusters for precision cosmology relies on accurate mea- 
surements of their gravitating mass and can be significantly af- 
fected by the presence of turbulence in the ICM. X-ray mass esti- 
mates are based on the assumption that the gas and the total mass 
distribution are spherically symmetric and in hydrostatic equi- 
librium (Rasia et al. 2006; Nagai, Vikhlinin & Kravtsov 2007; 
Jeltema et al. 2008; Piffaretti & Valdarnini 2008; Lau, Kravtsov 
& Nagai 2009). However, the presence of random gas motions 
implies additional pressure support which is not accounted for 
by the hydrostatic equilibrium equation. 

Turbulence in the ICM may also play an important role in 
non-thermal phenomena, such as the amplification of seed mag- 
netic fields via dynamo processes (Dolag, Bartelmann & Lesch 
2002; Subramanian, Shukurov & Haugen 2006) and the accel- 
eration of relativistic particles by magnetohydrodynamic waves 
(Brunetti & Lazarian 2007). The transport of metals in the ICM 
is also likely to be driven by turbulence (Rebusco et al. 2006). 

Numerical simulations provide a valuable tool with which to 
follow in a self-consistent manner the complex hydrodynamical 
flows which take place during the evolution of the ICM. In par- 
ticular, hydrodynamical simulations of merging clusters showed 
that moving substructures can generate turbulence in the ICM 
through shearing instabilities (Roettiger, Loken & Burns 1997; 
Norman & Bryan 1999a; Takizawa 2000; Ricker & Sarazin 
2001 ; Fujita, Suzuki & Wada 2004; Takizawa 2005; Dolag et al. 
2005; Iapichino & Niemeyer 2008; Vazza et al. 2009; Planelles 
&Quilis 2009). 

The solution of the hydrodynamic equations in a simulation 
depends on the adopted numerical method and in cosmology 
hydrodynamic codes which are used to perform simulations of 
structure formation can be classified into two main categories: 
Eulerian, or grid based codes, in which the fluid is evolved on 
a discretized mesh (Stone & Norman 1992; Ryu et al. 1993; 
Norman & Bryan 1999b; Fryxell et al. 2000; Teyssier 2002), 
and Lagrangian methods in which the fluid is tracked following 
the evolution of particles of fixed mass (Monaghan 2005). Both 
of these methods have been widely applied to investigate the for- 



mation and evolution of galaxy clusters (cf. Borgani & Kravtsov 
2009, and references therein). 

The main advantage of Lagrangian or smoothed particle hy- 
drodynamics (SPH) codes (Hernquist & Katz 1989; Springel, 
Yoshida & White 2001; Springel 2005; Wetzstein et al. 2009) 
is that they can naturally follow the development of matter con- 
centration, but they have the significant drawback that in order to 
properly model shock structure they require in the hydrodynamic 
equations the presence of an artificial viscosity term (Monaghan 
&Gingold 1983). 

Eulerian schemes such as the Parabolic Piecewise Method 
(Colella & Woodward 1984) implemented in ENZO (Norman 
& Bryan 1999b; O'Shea et al. 2005a) and FLASH codes 
(Fryxell et al. 2000), are instead characterized by the lack of 
artificial viscosity and by a better shock resolution when com- 
pared against SPH codes. The development of adaptative mesh 
refinement (AMR) methods, in which the spatial resolution of 
the Eulerian grid is locally refined according to some selection 
criterion (Berger & Colella 1989; Kravtsov, Klypin & Khokhlov 
1997; Norman 2005), led to a substantial increase in the dy- 
namic range of cosmological simulations of galaxy clusters and 
a better capability in following the production of turbulence in 
the ICM induced by merger events (Iapichino & Niemeyer 2008; 
Maier et al. 2009; Vazza et al. 2009; Vazza, Gheller & Brunetti 
2009). 

In principle, application of different codes to the same test 
problem with identical initial setups should lead to similar pre- 
dictions. However, a comparison between the results produced 
by AMR and SPH codes in a number of test cases reveals sev- 
eral differences (Frenk et al. 1999; O'Shea et al. 2005b; Agertz 
et al. 2007; Wadsley, Veeravalli & Couchman 2008; Tasker et 
al. 2008; Mitchell et al. 2009). Agertz et al. (2007) showed 
that the formation of fluid instabilities is artificially suppressed 
in SPH codes with respect AMR codes because of the difficul- 
ties of SPH codes in properly modeling the large density gradi- 
ents which develop at the fluid interfaces. The problem has been 
re-analyzed by Wadsley, Veeravalli & Couchman (2008), who 
concluded that the origin of the discrepancies is due partly to the 
artificial viscosity formulation implemented in SPH and mainly 
to the Lagrangian nature of SPH, which inhibits the mixing of 
thermal energy (Price 2008). 

Moreover, a long standing problem between the two numer- 
ical approaches occurs in non-radiative simulations of a galaxy 
cluster, where a discrepancy is produced at the cluster core by 
the two codes in the radial entropy profile of the gas. (Frenk et 
al. 1999; O'Shea etal. 2005b; Wadsley, Veeravalli & Couchman 
2008; Mitchell et al. 2009). To investigate the origin of this dis- 
crepancy, Mitchell et al. (2009) compared the final entropy pro- 
files extracted from simulation runs of idealized binary merger 
cluster simulations. They found that in the cluster central re- 
gions ( ~ few per cent of the virial radius) the entropy profile 
of Eulerian simulations is a factor ~ 2 higher than in the SPH 
runs. The authors argue that the main source of this difference in 
the amplitude of central entropy is strictly related to the amount 
of mixing present in the two codes. 

In the Eulerian codes it is the numerical scheme which forces 
the fluids to be mixed below the minimum cell size. This is in 
contrast with SPH simulations, in which some degree of fluid 
undermixing is present owing to the Lagrangian nature of the 
numerical method. When compared with the results of clus- 
ter simulations performed with AMR codes, entropy generation 
through fluid mixing is therefore inhibited in SPH simulations. 
Although a certain degree of overmixing could be present in 
mesh-based codes (Springel 2009), it appears then worth pur- 
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suing any improvement in the SPH method which leads to an 
increase in the amount of mixing present in SPH cluster simula- 
tions. This is motivated by the strong flexibility that Lagrangian 
methods posses in tracking large variation in the spatial extent 
and in the density of the simulated fluid. Fluid mixing is ex- 
pected to increase if viscous damping of random gas motions if 
effectively reduced. This effect occurs in SPH simulations be- 
cause of the presence of an artificial viscosity term in the hydro- 
dynamic equations, which is necessary to properly model shocks 
but introduces a numerical viscosity. 

In the standard formulation of SPH, the strength of artificial 
viscosity for approaching particles is controlled by a pair of pa- 
rameters which are fixed throughout the simulation domain. This 
renders the numerical scheme relatively viscous in regions far 
away from the shocks, with subsequent pre-shock entropy gen- 
eration and the damping of turbulent motions as a side effect. 
Given these difficulties Morris & Monaghan (1997) proposed a 
modification to the original scheme in which each particle has 
its own viscous coefficient, whose time evolution is governed by 
certain local conditions which depends on the shock strength. 
The benefit of this method is that the artificial viscosity is high 
in supersonic flows where it is effectively needed but quickly 
decays to a minimum value in the absence of shocks. As a con- 
sequence, the numerical viscosity is strongly reduced in regions 
away from shocks and the modeling of turbulence generated by 
shearing motions is greatly improved. 

Given the shortcomings of SPH codes which affect the de- 
velopment of turbulence in hydrodynamic cluster simulations, it 
is therefore interesting to conduct a study to analyze the effect 
that the adoption of a time-dependent artificial viscosity has on 
ICM properties of SPH cluster simulations. This is the goal of 
the present work, in which a time-dependent artificial viscosity 
formulation is implemented in a SPH code with the purpose of 
studying the differences induced in the ICM random gas motion 
of simulated clusters with respect the standard artificial viscosity 
scheme. 

To this end, a test suite of simulations is presented in which 
the ICM final profiles and turbulent statistical properties provide 
the quantitative measures used to compare runs with different 
artificial viscosity parameters, cluster dynamical states, numeri- 
cal resolution and physical modeling of the gas. A similar study 
has already been undertaken in pioneering work by Dolag et al. 
(2005), who analyzed the role of numerical viscosity and the 
level of random gas motions in a set of cosmological SPH clus- 
ter simulations. Here several aspects of a time-dependent arti- 
ficial viscosity implementation in SPH are discussed in a more 
systematic way. 

Specifically, we study the effects on the development of tur- 
bulence in the ICM which come from varying the simulation 
parameters which govern the time evolution of the artificial vis- 
cosity strength in the new formulation. Moreover, in addition to 
the adiabatic gas dynamical simulations a set of cooling runs was 
constructed in which the physical modeling of the gas includes 
radiative cooling, star formation and energy feedback. 

The paper is organized as follows. In Sect. 2 we present the 
hydrodynamical method and the implementation of the artificial 
viscosity scheme. The construction of the set of simulated cluster 
samples used to perform comparisons between results extracted 
from SPH runs with different artificial viscosity parameters, is 
described in Sect. 3. Sect. 4 provides an introduction to sev- 
eral statistical methods used in homogeneous isotropic turbu- 
lence to characterize statistical properties of the velocity field 
of a medium. Some numerical tests are discussed in Sect. 5, in 
order to assess the validity of the code and its shock resolution 



Table 1. Summary of the AV parameters used in the simulations. First 
column is the index label which identifies the runs with different AV 
parameters. The index is for the standard AV scheme, other indices 
refer to different values of the (a min , Id) pair used in the new time- 
varying AV scheme. The last column gives the value of decay constant 
S a ~ 0.447// d defined by Morris & Monaghan (1997). 



AV 






U 







l. 


l. 






1 


0.1 


1.5 


0.1 


4.5 


2 


0.1 


1.5 


0.2 


2.25 


3 


0.1 


1.5 


0.5 


1 


4 


0.1 


1.5 


1 


0.5 


5 


0.01 


1.5 


1 


0.5 



capabilities. The results of the cluster simulations are presented 
in Sect. 6, while Sect. 7 summarizes the main conclusions. 



2. Description of the code 

Here we provide the basic features of the numerical scheme used 
to follow the hydrodynamics of the fluid; for a comprehensive 
review see Monaghan (2005). 



2. 1 . The hydrodynamical method 

The hydrodynamic equations of fluid motion are solved accord- 
ing to the SPH method, in which the fluid is described within 
the domain by a collection of N particles with mass m„ veloc- 
ity V,-, density p, and a thermodynamic variable such as the spe- 
cific thermal energy m, or the entropy A,. The latter is related 
to the particle pressure P, by P, = A,pJ, where y — 5/3 for a 
monoatomic gas. The density estimate p(r) at the particle posi- 
tion r, is given by 



Pi = ^m j W(r i j,h i ) , 



(1) 



where W(|r,- - rj\,h,) is the B2 or cubic spline kernel which has 
compact support and is zero for |r,-r,| > 2h, (Monaghan 2005). 
The sum (1) is over a finite number of particles and the smooth- 
ing length hi is a varying quantity which is implicitly defined 
through the equation (Springel & Hernquist 2002) 



47r(2/z,) 3 Pi 



: Nsphfiii , 



(2) 



with N sp h being the number of neighboring particles within a 
radius 2hi. Typical choices lie in the range N sp h ~ 33 - 50, here 
N sp h = 33 is used. 

The equation of motion for the SPH particles is given by 



dVi vh 



Pi 



V l W u (h i ) + fj-+ 1 V i W ii (hj) 



where the coefficients £1 are defined as 



£1 = 



Shi sr dW ik (hd 



Olli ST- 



(3) 



(4) 



and in the momentum equation (3) account for the effects due to 
the gradients of the smoothing length hi (Monaghan 2005). The 
momentum equation must be generalized by including an addi- 
tional viscous pressure term which in SPH is needed to represent 
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Table 2. Main cluster properties and simulation parameters of the baseline sample. From left to the right : the dynamical subset class , the 
cosmological sample membership, the corresponding sample index of the cluster, the present cluster mass M 2 oo within r 2 oo in units of hT l M Q , the 
radius r 2 oo is units of Mpc, the number of gas particles N gas inside the L c /2 sphere at z = Zi„, the same but for the dark matter particles, the mass of 
the gas particles m gas in M Q , the same but for the dark matter particles, the gravitational softening parameter for the gas in kpc and the root mean 
square plane average cluster power ratio n 3 (r 20 o) 
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the effects of shocks. This is achieved in SPH by introducing an 
artificial viscosity (herefater AV) term with the purpose of con- 
verting kinetic energy into heat and preventing particle interpen- 
etration during shocks. The new term is given by 



'dVi 
t dt 



/ vise ; 



(5) 



where the term Wij = |(W(r,y, hi) + W(rij,hf)) is the sym- 
metrized kernel and n, 7 is the AV tensor. To follow the thermal 
evolution of the gas, an entropy-conserving approach (Springel 
& Hernquist 2002) is used here and entropy is generated at a 
rate 



PI 



(6) 



where Vy 



v, - Vj, Ti is the particle temperature and the 
additional term A(p,, Tf) accounts for the radiative losses of the 
gas, if present. In the following, simulations in which the cooling 
term A is absent from Eq. (6) will be referred to as adiabatic. The 
expression for the AV tensor If,-, is 



-aijdjUij+PijH^ 



n, 







Pij 



fa if v u ■ r u > 
otherwise , 



(7) 



and so n, ; - is non-zero only for approaching particles. Here scalar 
quantities with the subscripts ; and j denote arithmetic averages, 
d is the sound speed of particle i, the parameters a, and reg- 
ulate the amount of AV and is a controlling factor which re- 
duces the strength of AV in presence of shear flows. The f^j term 
is given by 



W; = 2 



hifij ■ r u 



(8) 



where the factor 77 = 10~ 2 is included to prevent numerical di- 
vergences. In order to limit the amount of AV generated in shear 
flows, Balsara (1995) proposed the following expression for : 



fi = 



IV • v\i 



V • v\i + |V x v\i + rf. ' 



(9) 



where (V • v), and (V x v), are the standard SPH estimates for 
divergence and curl (Monaghan 2005), and the factor v[ i = 
10~ 4 Ci/hj is inserted to prevent numerical divergences. This ex- 
pression for is effective in suppressing AV in pure shear flows, 
for which the condition |V x v|, » |V • v|, holds. 



The strength of the AV in the standard SPH formulation is 
given by = 2a i and or,- = const = ao, with ao = 1 being a 
common choice for the viscosity coefficient (Monaghan 2005). 
In the following this parametrization will be referred to as the 
'standard' AV. 

Monaghan (1997) proposed a modification to this 
parametrization for AV based on an analogy with the Riemann 
problem. In a number of test problems, results obtained us- 
ing his 'signal velocity' formulation are found to be equivalent 
or slightly improved with respect to the standard AV scheme. 
However, this new formulation is not introduced here in order 
to avoid any further source of difference in the results produced 
by the SPH code, additional to those due to the choice of the 
time-dependent AV scheme parameters. 

In simulations in which the gas is allowed to cool radiatively, 
cold gas in high density regions is subject to star formation and 
gas particles are eligible to form star particles. Energy and metal 
feedback is returned from these star particles to their gas neigh- 
bors by supernova explosions, according to the stellar lifetime 
and initial mass function. For a detailed description see also 
Valdarnini (2006). 

2.2. The new artificial viscosity formulation 

The standard AV formulation is successfull in properly resolv- 
ing shocks but at the same time generates unwanted viscous 
dissipation in regions of the flow which are not undergoing 
shocks. Morris & Monaghan (1997) proposed that the viscous 
coefficient should be different for each particle and left free to 
evolve in time under the local conditions. Following Morris & 
Monaghan (1997) 



dot on - a„ 



dt Ti 
where 

hi 

Ci k 



(10) 



(11) 



is a decay time-scale which regulates, through the dimensionless 
parameters I4, the time evolution of a,(f) away from shocks. The 
parameter a min is the minimum value to which a,(?) is allowed 
to decay. The source term S , is given by 

Si = fiS max(-(V • v),-, 0)(a max - a t ) = S i(a max - a t ) , (12) 

and is constructed in such a way that it increases in the pres- 
ence of shocks. Following a suggestion of Morris & Monaghan 
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(1997), the damping factor f t is inserted to account for the pres- 
ence of vorticity. The scale factor 

5/3+1 y+1 
5 = ln(^— r)/ln(— r) (13) 
5/3 - 1 y - 1 

is chosen such that the peak value of a, is independent of the 
equation of state. The source term equation (12) is of the form 
proposed by Rosswog et al. (2000), which allows a better re- 
sponse in the presence of shocks with respect to the original for- 
mulation. In a number of test simulations Rosswog et al. (2000) 
found that appropriate values for the parameters a max , a min , Id are 
1.5, 0.05 and 0.2, respectively. 

From the viewpoint of the description of the fluid flow veloc- 
ities, the most significant parameters are a m ,„ and Id- Since the 
goal of the SPH simulations presented here is to investigate the 
effect on ICM fluid flows of the numerical viscosity, the set of 
runs will be performed with the parameter a max set to the value 
<Xmax =1.5 (see Sect. 5.1.1), whilst a range of values will be con- 
sidered for the AV parameters a m ,„ and Id- However, a lower limit 
to the timescale t,- is set by the minimum time taken to propagate 
through the resolution length /z„ so that the value Id = 1 sets an 
upper limit to the parameter Id- 

The different implementations of AV used in the simulations 
are summarized in Table 1. In the simulations which incorpo- 
rate the new time-dependent AV scheme, five different pairs of 
values have been chosen for the parameters (a m! „, Id), while sim- 
ulations in which the AV is modeled according to the standard 
formulation are used for reference purposes. Simulations with 
different AV schemes or parameters will then be labeled by the 
corresponding index of Table 1 . 

However, if the decay parameter Id approaches unity, the 
time-dependent AV scheme discussed here may fail to prop- 
erly evolve the viscosity parameters a, when a shock is present. 
To avoid these difficulties an AV scheme which generalizes the 
Rosswog et al. (2000) source term expression is proposed here. 
Since the implementation of the new source term equation is 
strictly related to the shock tube problem discussed in Sect. 5.1, 
the modifications to Eq. (12) will be presented in Sect. 5.1 to- 
gether with the shock tube tests. 

3. Sample construction of simulated clusters 

In order to investigate the effects of numerical viscosity on ICM 
random gas motions a large ensemble of hydrodynamical clus- 
ter simulations was created by performing, for a chosen AV test 
case, SPH hydrodynamical simulations using a baseline sample 
consisting of eight different initial conditions for the simulated 
clusters. These initial conditions are those of simulated clusters 
which are part of a large set produced in an ensemble of cosmo- 
logical simulations (see later), and these clusters are chosen at 
the present epoch with the selection criterion of covering a wide 
range in virial masses and dynamical properties. This choice was 
made in order to study, using the hydrodynamical simulations, 
the dependencies of the ICM gas velocities on these cluster prop- 
erties. The number of eight clusters was a compromise between 
the need for obtaining results of sufficient statistical significance 
and of keeping the computational cost at a minimum level, given 
the range of parameters explored. 

The initial conditions of the baseline sample were con- 
structed according to the following procedures (Valdarnini 
2006), a similar approach has been followed by Dolag et al. 
(2005) to construct their set of high resolution hydrodynami- 
cal cluster simulations. An N-body cosmological simulation is 



first run with a comoving box of size La = 200/i~ 1 Mpc. The 
cosmological model assumes a flat geometry with the present 
matter density Q m = 0.3, vacuum energy density Q A = 0.7, 
£l h = 0.0486, and with h = 0.7 being the value of the Hubble 
constant H in units of 100 km s~'Mpc~'. The scale invariant 
power spectrum is normalized to cr 8 = 0.9 on an 8 hr l Mpc scale 
at the present epoch. 

A cluster catalog was generated by identifying clusters in the 
simulation at z — using a friends-of-friends (FoF) algorithm to 
detect overdensities in excess of ~ 200O m ° 6 . The sample com- 
prises N2 = 120 clusters ordered according to the value of their 
mass M200, where 

M A = (4n/3)A Pc rl (14) 

denotes the mass contained in a sphere of radius r& with 
mean density A times the critical density p c . 

Each of these clusters was then resimulated individually us- 
ing an N-body+SPH simulation performed in physical coordi- 
nates starting from the initial redshift z,„ = 49. The integration 
was made by first locating the cluster center at z = and iden- 
tifying all of the simulation particles lying within r2oo- These 
particles are then located at z,„ and a cube of size L c oc M^l , 
enclosing all of them, is found. A lattice of Nl = 74 3 is set in- 
side the cube and to each node is associated a gas particle with 
its mass and position, a similar lattice is set for dark matter par- 
ticles. The particle positions are then perturbed, using the same 
random realization as for the cosmological simulations. The par- 
ticles kept for the hydrodynamic simulation are those whose per- 
turbed positions lie within a sphere of radius L c /2 from the clus- 
ter center. To model tidal forces, the sphere is surrounded out 
to a radius L c by an external shell of dark matter particles ex- 
tracted from a cube of size 2L C centered in the same way as 
the original cube and consisting of Nl = 74 3 grid points. The 
particles are evolved up to the present time using an entropy- 
conserving SPH code, described in Sect. 2.1, combined with a 
treecode gravity solver. Particles are allowed to have individual 
timesteps and their gravitational softening parameter is set ac- 
cording to the scaling e, am/ , where m, is the mass of particle 
;. The relation is normalized by e, = 15 (m,/6.2 • 10 8 M o )^ 3 kpc. 
The set of these individual cluster hydrodynamical simulations 
is then denoted as sample S2- 

The whole procedure is then repeated twice more in order to 
generate the cluster samples S 4 and from cosmological sim- 
ulations with box size L 4 = 400/r'Mpc and Lg = 800/z _1 Mpc. 
The number of clusters N m in these samples is chosen such that 
the mass M 2 oo of the N,„ - th cluster of sample S m is greater 
that the mass M200 of the most massive cluster of sample S m /2, 
and samples S& and S4 consist of N& = 10 and N4 = 33 clus- 
ters, respectively. The final cluster sample S a ii is constructed by 
combining all of the samples S m generated from the three cos- 
mological runs. 

In order to provide the baseline sample for the hydrodynami- 
cal simulations, eight clusters were chosen from those of sample 
S a u, with the selection criterion being the construction of a rep- 
resentative sample of the cluster dynamical states and masses. 

As an indicator of the cluster dynamical state, we adopt the 
power ratios P r /Po- The quantity P r is proportional to the square 
of the r-th moments of the X-ray surface brightness S x(x,y), as 
measured within a circular aperture of radius R ap in the plane 
orthogonal to the line of sight (Buote & Tsai 1995). Here the 
power ratio P3/P0, or equivalently rii,(R d p) = log lQ (Pj,/Po), is 
used as an estimator of the cluster dynamical state since it gives 
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Table 3. The main parameters used in simulations of different resolu- 
tion. The first column gives the label which denotes the runs with dif- 
ferent numbers of particles, N L is the number of grid points set inside 
the L c cube, and the last two columns indicate the approximate number 
of gas particles and total particles used in the cluster simulations. 
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an unambiguous detection of asymmetric structure. For a fully 
relaxed configuration Il r — > -oo. 

For the simulated clusters of sample S a u, the power ratios 
are evaluated at z = in correspondence of the aperture radius 
Rap = ?"200- To minimize projection effects, the average quantity 
Fb^oo) = log 10 (P3/Po) is used to estimate the cluster dynami- 
cal state, where P r it is the rms plane average of the moments P r 
evaluated along the three orthogonal lines of sight. Clusters of 
sample S a ii are then sorted according to their values of Fbfooo) 
and eight clusters are extracted from the sample. The four of 
them with the lowest values of FFt^oo) among the cluster power 
ratio distribution are chosen and these relaxed clusters will be 
denoted as quiescent (Q) or relaxed clusters. The remaining four 
are chosen with the opposite criterion of having the power ra- 
tios with the highest values among those of the 03(^200) distri- 
bution. These clusters are then labeled as perturbed or P clusters. 
Within each subset, clusters are chosen with the additional crite- 
rion of having their virial mass distribution as wide as possible. 
The cluster mass at r2oo is used in place of the virial mass and 
Table 2 lists the main cluster properties and simulation param- 
eters of the baseline sample constructed according to the above 
criteria. 

The numerical convergence of the results is assessed by 
studying the dependence of the final profiles on the numerical 
resolution adopted in the simulations. In addition to the hydro- 
dynamical simulations realized from the baseline sample using 
various AV prescriptions, a set of mirror runs with a different 
number of simulation particles was then performed with the aim 
of studying the stability of the final results against the numerical 
resolution of the simulations. Because of the large number of AV 
test simulations performed here, the stability against resolution 
of the simulation results for the baseline sample is investigated 
by comparing the chosen profiles with the corresponding ones 
obtained from parent simulations with lower resolution. A suf- 
ficient condition for the numerical convergence of the profiles 
is given by their stability against simulations performed with 
higher resolution. However, the approach undertaken here pro- 
vides significant indications about the stability of the baseline 
sample and at the same time allows a substantial reduction in the 
amount of computational resources needed to assess numerical 
convergence. 

The initial conditions of the mirror runs are therefore con- 
structed using the same procedures as previously described for 
the baseline sample, the only difference being the number Nl 
of grid points set inside the L c cubes. With respect to the base- 
line sample, which will be denoted as high resolution (HR), low 
(LR) and medium resolution (MR) runs are obtained by setting 
Nl = 35 3 and Nl = 51 3 , respectively. Table 3 lists some basic 
parameters. Individual runs have their label obtained by combin- 
ing the different labels specified in the Tables. 



4. Statistical measures 

Here we present several statistical tools with which will be stud- 
ied the impact that the strength of AV has on the statistical prop- 
erties of the turbulent velocity field u(x) of the simulated clus- 
ters. Homogeneous isotropic turbulence is commonly studied us- 
ing the spectral properties of the velocity field u(x). Its Fourier 
transform is defined as 



(2tt) j J 



(15) 



and the ensemble average velocity power spectrum V(k) is given 
by 



< u l (k')u(k) >= 5 D (k' - k)P(k). 

The energy spectrum function Eik) is then defined as 

E(k) = 2nk 2r P(k), 



(16) 



(17) 



where k = \k\, and the integral of E(k) is the mean kinetic energy 
per unit mass 



(18) 



E(k)d k= - <u 2 > . 
o 2 

In the case of incompressible turbulence the energy spec- 
trum follows the Kolgomorov scaling E(k) oc k~ 5/3 . In contrast, 
the regime of compressible turbulence exhibits large variation in 
the gas density and a generalized energy spectrum can be con- 
sidered by introducing in Eq. (15) a weighting function u(x) — > 
u w {x) = w(x)u(x), where w(x) is proportional to some power 
of the density. For compressible turbulence a natural choice is 
w(x) oc yfp(xj, in this case the integral of E(k) is just the kinetic 
energy density : J E(k)d k = j < pu 2 >. As noticed by Kitsionas 
et al. (2009) the spectral properties of compressible turbulence 
are better described using this weighting scheme, which allows a 
more accurate estimate of the power at small scales where mass 
accumulation occurs because of shocks. Alternatively, Kritsuk 
et al. (2007) demonstrated that for supersonic isothermal turbu- 
lence, a Kolgomorov scaling for the spectrum function E(k) is 
recovered using w(x) oc p(x) 1 ^ 3 . Here the spectral properties of 
compressible turbulence will be studied using a generalized en- 
ergy spectrum with density weighting given by w(x) oc p(x) 1/2 . 
This is the same choice adopted by Kitsionas et al. (2009) and 
provides for the energy spectrum E(k) a physical reference to the 
kinetic energy density. 

To study the energy spectrum in cluster simulations, spec- 
tral quantities were constructed as follows. A cube of side L sp 
with Ng grid points is first placed at the cluster center and hy- 
drodynamic variables A(x) are estimated at the grid points x p 
according to the SPH prescription 



A(x p )=Y j A i ^W(\x p -x i \,h i ), 



(19) 



where A(x) denotes either the components of the velocity 
field u(x) or the weighting function w(x). 

The weighting function is set to w(x) = p(x) l/2 /S w , where 
S w = J^p(x p ) l ^ 2 /Ng. This normalization choice guarantees 
2 w(x p ) = Ng, regardless of the power density exponent used 
in the weighting. The discrete transforms of u„(x p ) are then 
computed using fast Fourier transforms and an angle-averaged 
density-weighted power spectrum f^QC) =< \u w d (k)\ 2 > is gen- 
erated as a function of k — \k\, binning the quantity \u w d {k)\ 2 in 
spherical shells of radius k and averaging in the bins. 
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Fig. 1. Results at t = 0. 12 from the shock tube test for 3D SPH runs with different AV parameters. The profiles are projected along the shock front. 
From top to bottom are plotted: density, velocity and thermal energy ( left column); pressure, entropy and viscosity parameter (right column). The 
solid black line represents the analytical solution, while lines with different styles and colors are the profiles of the SPH runs with different AV 
parameters. Different runs are labeled according to Table 1 and the relationship with the corresponding profiles is illustrated in the entropy panel. 
In the density panel profiles of different runs have been shifted vertically to better illustrate their relative differences. 



As an estimator of the the energy density (17) one can there- 
fore use E(k) = 2jik 2( P d {k)(^-) i . However, in order to consis- 
tently compare density-weighted power spectra extracted from 



different clusters and boxes it is useful to define a rescaled tur- 
bulent power spectrum as 



E(k) = 



1 



Lsncrl 



In 



(20) 
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where cr v = VGA^oo/^oo- The dimensionless form of this 
power spectrum allows one to compare curves of E(k) extracted 
from different clusters as a function of k = k r L sp /2n, where 
k r = \k\. The spectral properties of the gas velocity fields will 
be studied using the density-weighted power spectrum defined 
according to Eq. (20), whereas volume-weighted (w — 1) spec- 
tra will be considered for comparative purposes. 

Moreover, the longitudinal and solenoidal components of the 
power spectrum E(k) will be studied separately. For doing this, 
the shear and compressive parts of the velocity u{x) are defined 
respectively in the £-space as 

k x u(k) 

u{k) shear = — — — , (21) 
1*1 

k ■ u(k) 

M(.fc)comp ~ jT"j ■ (22) 

The corresponding power spectra are used to define the 
power spectrum decomposition, E(k) = E s (k) + E c {k) (Kitsionas 
et al. 2009; Zhu, Feng & Fang 2010). 

In addition to the spectral properties of the velocity u(x), it is 
useful to investigate the scaling behavior of the velocity structure 
functions. These are defined by 

S p (r) =< \u(x + f) - u(x)\ p >, (23) 

where p is the order of the function. For computing the structure 
functions, a random subsample of N s particles is constructed; 
these particles are extracted from the set of gas particles which 
satisfy the constraint of being located within a distance r^oo from 
the cluster center. For each of these sample particles s, the rel- 
ative velocity difference Ah = u(Xj + r s i) - u(x s ) is computed 
for all of the gas particles i separated by a distance r s , < r2oo, 
and the quantity |Am| p is binned in the corresponding radial bin. 
Final averages are obtained by dividing the binned quantities by 
the corresponding number of pairs belonging to the radial bin. 

As for the power spectrum, structure functions have been 
computed separately for both transverse and longitudinal com- 
ponents. These are respectively defined as Au ± = Am x r s ;/|r s ;| 
and A«|| = Am • r S i/\r s i\. Here the study will be restricted to 
second-order (p = 2) structure functions. Moreover, in order 
to compare structure functions of different clusters, these will be 
rescaled according to S p (r) - S p (r)/a^ and radial distances will 
be expressed in units of r2oo. Density- weighted velocity structure 
functions are defined using the weighted field u w (x) in Eq. (23). 

Finally, the energy content of the turbulent velocity field is 
another quantity used to investigate the dependence of the veloc- 
ity field properties on the amount of AV present in the simula- 
tions. In order to properly consider the contribution of random 
gas motion to the turbulent energy budget, it is however neces- 
sary to separate from the velocity field u(x) the part due to the 
streaming motion. A useful approach is then to define a spatial 
decomposition which separates the velocity into a small-scale 
and a large-scale component (Adrian, Christensen & Liu 2000). 
More specifically, 

<u(x,f)>= f f(x-x',H)u(x',f)d 3 x' , (24) 
Jd 

where f(x - x',H) is a shift-invariant kernel which acts as 
a low-pass filter, H being a filtering scale. Velocity fields are 
always considered here at a given time slice, so that the depen- 
dence on t will be removed from Eq. (24). The small-scale tur- 
bulent velocity field u{x) is then defined as 

u(x) = u(x)- < u(x) > . (25) 
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Fig. 2. A closer view at the shock front of the density profiles displayed 
in the top left panel of Fig. 1. 

This method of decomposing a turbulent velocity field is also 
known as Reynolds decomposition and is commonly employed 
in Large Eddy Simulations (LES, Meneveau & Katz 2000). 
The filtering kernel chosen here to define the large-scale veloc- 
ity field is a triangular-shaped cloud function (TSC) (Hockney 
& Eastwood 1988). This choice is similar to that adopted by 
Dolag et al. (2005) when studying the ICM velocity fields, and 
is motivated by the compact support properties of the kernel, 
the TSC being a second-order scheme, which from a computa- 
tional viewpoint greatly reduces the complexity of estimating the 
mean field velocity (24). The choice of the filtering scale H for 
the TSC spline will be discussed in Sect. 6; here is worth noting 
that after application of Eq. (24) to u(x), the harmonic content 
of the small-scale turbulent velocity field u{x) is defined by the 
condition kH » 1. 

5. Test simulations 

In this section simulation results are presented obtained by ap- 
plying the SPH implementation described in the previous sec- 
tions to several test problems. These tests have known solutions 
provided by analytical models or independent numerical meth- 
ods, and against these solutions are compared results from the 
SPH runs with the purpose of validating the code and assessing 
its performance. Section 5.1 is dedicated to the shock tube prob- 
lem and Sect. 5.2 to the 3D collapse of a cold gas sphere. Both 
of these tests have been chosen because they have been widely 
used in the literature and moreover they allow comparisons with 
previous SPH runs in which a time-dependent AV scheme was 
implemented. 

5.1. The shock tube problem 

5.1 .1 . Initial condition set-up and parameter calibration of the 
time-dependent AV scheme 

The Riemann shock tube problem is a test commonly used for 
SPH codes. Its main advantage is that it admits an analytical so- 
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lution following the propagation of a shock wave in a medium 
initially at rest. Here the initial condition setup consists of an 
ideal fluid with y = 5/3 initially at rest at t = 0. An in- 
terface at x = separates the fluid on the right with density 
and pressure (pi,P\) = (4, 1) from the fluid on the left with 
(P2, Pi) = (1, 0.1795). The fluid is then allowed to evolve freely 
and at later times there is a shock wave propagating toward the 
left, followed by a rarefaction wave and a contact discontinuity. 
The resulting profiles are given by the analytic solution (Rasio 
& Shapiro 1991). 

The one-dimensional Riemann shock tube problem is often 
used to test hydrodynamic codes, but the numerical results are 
often of limited validity because numerical effects which may 
arise in 3D calculations, such as particle streaming, are absent 
or reduced when there is only one degree of freedom. For this 
reason the shock test is carried out here in three dimensions. To 
construct the initial conditions for the SPH runs, a cubic box of 
side-length unity was filled with 10 6 equal mass particles. Of 
these one million particles, 800 000 were placed in the right-half 
of the cube, while 200 000 were placed in the left-half. The parti- 
cles in the two halves of the cube were extracted from two inde- 
pendent uniform glass-like distributions contained in a unit box. 
The two distributions had 1.6 • 10 6 and 4 • 10 5 particles, respec- 
tively. 

These initial conditions have been chosen because they are 
the same as those adopted by Tasker et al. (2008, hereafter T08), 
in their Sect. 3. 1, to study the shock tube problem. Those authors 
used a variety of numerical problems with known analytical so- 
lutions to compare the behavior of different astrophysical codes 
and their ability to resolve shocks. It is therefore of particular 
interest to compare the results of T08 with those obtained from 
the shock tube SPH simulations performed here using different 
AV strengths. The SPH runs were realized imposing periodic 
boundary conditions along the y and z axes and the results were 
examined at t = 0.12. Runs with the time-dependent AV scheme 
have their viscosity parameters initialized to one. 

Before proceeding to discuss the behavior of the shock tube 
tests, it is necessary to assess the validity, under different condi- 
tions, of the AV scheme introduced in Sect. 2.2. As outlined in 
this section, the time evolution of the viscosity parameter can be 
affected if very short damping time scales are imposed. To bet- 
ter illustrate this point it is useful to rewrite Eq. (10) using the 
second of the equalities of Eq. (12). The new Eq. (10) is then 



dt 

where 



T; = 



Ti 



and qt is a modified source term 



l+Sm 



(26) 



(27) 



(28) 



Neglecting variations in the coefficients, the solution to Eq. 
(26) at times t > t m is given by 



«r/(0 = qi + {aAtin) - qd exp {t '"' )lT '> 



(29) 



which shows that a,(f) approaches the asymptotic value a m! „ in 
the absence of shocks, and a max if a shock is present. However, 
for the source term qi, the condition qi a max holds only in 
the shock regime S,t, » 1. If this condition is not satisfied and 



SiTi £ 1, the peak value a pea k of a,(/) will be smaller than a max . 
In particular a pea k will depend on the chosen value of the decay 
parameter la- This dependence of a pea k on la is absent for strong 
shocks, but for mild shocks introduces the unwanted feature that 
for short decay timescales (la — > 1) the peak value of ar,-(f) at 
the shock front might be below the AV strength necessary to 
properly treat shocks. 

In fact, application of the SPH code implemented according 
to the AV scheme described in Sect. 2.2 to the shock tube prob- 
lem with the initial condition setup studied here, shows that at 
the shock front location the peak in the AV parameter ranges 
from unity down to ~ 0.2 for l d = \. Previous numerical tests 
(Morris & Monaghan 1997; Rosswog et al. 2000) showed that 
for this shock tube problem there is good agreement with the ex- 
act results using la = 0.2, and at the shock front the peak in the 
viscosity parameter is then approximately ~ 0.6 - 0.7. 

In order to mantain the same shock resolution capabilities 
also in those cases studied using short time scales with la = 1, 
a correction factor £ is introduced so as to compensate in the 
source term q, the smaller values of S ,t, with respect to the small 
Id regimes. This is equivalent to considering a higher value for 
a max , so that in Eq. (28) a max is substituted by a max -> £a max . 
The calibration of the correction factor £ is achieved by requir- 
ing that for the current shock tube problem, the SPH runs with 
AV decay parameters l d > 0.2 should have at the shock front the 
same value a pea k as that for the la = 0.2 case, i.e. a pea k ~ 0.6-0.7 
for la > 0.2. This choice guarantees that the viscosity parame- 
ter reaches the correct size at the shock front, but away from it 
quickly decays according to the chosen value of Id- After several 
tests it has been found that the choice £ = MAX((l d /0.2)° 8 , 1) 
yields satisfactory results and all of the numerical tests shown in 
this paper incorporate in the source term (12) the modification 



i = MAX((l d IQ.2f\ 1). 



(30) 



Note that the validity of this setting has been established for 
the shock strength of the shock tube problem considered here, 
nonetheless the results of the other tests indicate that this 
parametrization is appropriate to make the peak value of the vis- 
cosity parameter at the shock location independent of the chosen 
value for the decay parameter la- 

5.1.2. Results 

The results of the shock tube tests are shown in Figures 1 and 
2. In each panel, different profiles refer to runs with different 
AV parameters. The solid line is the analytical solution, which 
exhibits a shock front at x = -0.095 and a contact discontinuity 
at x = -0.033. The profiles shown in Figures 1 and 2 can be 
compared with the corresponding Figures 2 and 3 of T08. 

The fixed AV simulation is in good agreement with the 
GADGET2 run of T08, although a closer view of the density 
profile in the proximity of the shock front in Fig. 2 shows post- 
shock ringing features which are absent in T08. This happens 
because close to the shock the amount of AV generated by the 
code allows particle interpenetration. Note that in order to avoid 
post-shock oscillations in the density, T08 explicitly made the 
choice of setting the viscosity parameter in the GADGET2 run 
to ArtBulkVisc=2, whereas here the fixed AV run was performed 
with the viscosity parameter a set to unity. The post-shock os- 
cillations in velocity are very similar in amplitude and behavior 
to the ones produced by the GADGET2 run of T08, and simi- 
larly for the glitch in pressure at the contact discontinuity. The 
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Fig. 3. Results at / = 0.9 for the adiabatic collapse of a cold gas sphere. Clockwise from the top left panel: radial density profiles of density, 
pressure, viscosity parameter a and radial velocity. The black solid lines are the results of the ID PPM calculation of Steinmetz & Mueller (1993). 
SPH runs with different AV parameters are labeled in the same way as in the shock tube test panels of Fig .1. 



spike in thermal energy and entropy at the contact discontinuity 
is originated by the initial discontinuity in the density profile, 
which has been left unsmoothed at the beginning of the simu- 
lation. These features are present also in T08, although here the 
height of the spikes is a bit smaller. 

The profiles of the simulations in which a time-dependent 
AV scheme has been implemented show a behavior very simi- 
lar to those of the fixed AV run. The post-shock oscillations in 
the density have a tendency to be amplified as the AV scheme 
uses shorter decay time scales, but the effect is minimal. This 
can be seen in the top left panel of Fig. 1, in which the density 
profiles have been shifted vertically to illustrate the effect better. 
The differences in the other profiles as a function of the AV co- 
efficients are very small, with the post-shock entropy being the 
quantity with the strongest dependence and the smallest values 
being those of the AV runs with the lowest viscosity parameter 
(Id = 1). Finally, the radial profiles of the time-dependent arti- 
ficial viscosity parameters a are consistent with the discussion 
of Sect. 2.2: their peaks are located approximately at the shock 
front and their post-shock radial decay is faster for those AV runs 
with the shortest decaying time scales. 



To summarize, the simulations of the shock tube problem 
performed with the SPH code presented here using the standard 
AV scheme show results in good agreement with the analytic so- 
lution and with the ones obtained by other authors using similar 
initial condition setup and parameters. The profiles of the runs 
with a time-dependent AV scheme agree well with those of the 
fixed AV run, thus indicating shock resolution properties com- 
parable to those of the standard AV scheme, but with a reduced 
AV strength elsewhere. 

5.2. Collapse of a cold gas sphere 

The Riemann shock tube test of the previous section illustrates 
the ability of the code to resolve discontinuities, but the shock 
strengths produced by the test are much smaller than those which 
develop during the gravitational collapse of astrophysical ob- 
jects. A more stringent test for SPH codes in which gasdynamics 
is modeled including self-gravity is therefore to study a 3D col- 
lapse problem. The test used by Evrard (1988) follows in time 
the adiabatic collapse of a gas sphere. This test has been widely 
used by many authors (Hernquist & Katz 1989; Steinmetz & 
Mueller 1993; Dave, Dubinski & Hernquist 1997; Wadsley, 
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Fig. 4. The same as in Fig. 3, radial profiles of density and entropy are shown in the proximity of the shock front. 



Stadel & Quinn 2004; Springel 2005; Vanaverbeke et al. 2009; 
Wetzstein et al. 2009) as one of the standard tests for SPH codes. 
The shock strength which develops during the collapse is compa- 
rable to that of the Navarro & White (1993) self-similar collapse 
test. 

The gas cloud is spherically symmetric and initially at rest, 
with mass M, radius R and density profile 



The gas obeys the ideal gas equation of state with y = 5/3 
and the thermal energy per unit mass is initially set to u — 
0.05GM/R. The chosen time unit is the cloud free-fall time scale 
t ff = (7r 2 /8) y/R 3 /GM 1 and the SPH simulations are performed 
using units for which G — M — R — 1. The cloud begins to 
collapse under its own gravity and the gas in the core is heated 
until the growth of pressure is sufficient for the infalling mate- 
rial to bounce back and a shock wave then propagates outwards. 
Most of the kinetic energy is converted into heat at the epoch of 
maximum compression of the gas, which occurs at t ~ 1.1. 

The initial conditions for the simulations are constructed by 
stretching the radial coordinates of a glass-like uniform distri- 
bution of N - 88 000 equal mass particles contained within a 
sphere of unit radius. Radial coordinates are transformed accord- 
ing to the rule r — > r' — r 3 ^ 2 , so as to generate the density profile 
(31). The particle smoothing lengths are adjusted according to 
Eq. (2) with N sp h = 50 neighbors and the gravitational softening 
length is taken as s g = 0.02. This choice for the simulation pa- 
rameters allows one to compare the test calculations performed 
here with the results presented for the same collapse problem by 
Wetzstein et al. (2009) in their Sect. 7.1. 

Fig. 3 shows the radial dependence at t — 0.9 of the spher- 
ically averaged profiles for density, pressure, radial velocity 
and time-dependent viscosity. The solid line is for the high- 
resolution ID PPM simulation of Steinmetz & Mueller (1993), 

1 Note that this time normalization is the same as that of Hernquist 
& Katz (1989) and differs by a factor from that of Evrard 

(1988) and by jt/S with respect to that adopted by Steinmetz & Mueller 
(1993). 



which for practical purposes can be considered as being an ex- 
act solution. The differences between the profiles with different 
AV parameters are minimal and the post-shock radial decay of 
the viscosity parameters is in accord with the expectation of the 
model. The profiles of the SPH runs are consistent with previous 
findings (Springel 2005; Wetzstein et al. 2009) and reproduce 
the overall features of the PPM solution, with some differences 
at the shock front, which at this epoch is located at r ~ 0.18. 
To illustrate these differences, Fig. 4 shows a closer view of the 
density and entropy profiles in the proximity of the shock front. 

A feature common to all of the runs is a significant amount of 
pre-shock entropy generation. This is inherent in the AV imple- 
mentation of the SPH code, which near to the bounce is switched 
on by the strong convergence of the flow, therefore generating 
dissipation (Wadsley, Stadel & Quinn 2004). Clearly, different 
settings of the AV parameters are reflected in the code's abil- 
ity to model shocks and compressions which in turn generates 
the corresponding differences in the pre-shock entropy profiles. 
All of the simulations share the basic attribute of having a shock 
front located outside the position indicated by the reference PPM 
entropy profile. The location of the shock front has a weak de- 
pendence on the implemented AV settings, the time-dependent 
AV runs with the shortest decay time scales (runs AV4 and AV5) 
being the closest to the PPM reference position. Moreover, the 
fixed AV run overpredicts the entropy peak whilst the other time- 
dependent AV simulations are in rough agreement with the PPM 
one. Finally, the post-shock entropy profile of all the runs is be- 
low the corresponding reference PPM profile. 

The density and entropy panels of Fig. 4 can be compared 
with Figures 7 and 8 of Wetzstein et al. (2009). For the chosen 
parameter settings, the AVj run corresponds to the Wetzstein et 
al. (2009) model with 6 a = 5 and = 0.1, whilst AV2 corre- 
sponds to the 5 a = 2, a+ =0.1 model. A comparison shows that 
the results presented here share some common properties with 
those of Wetzstein et al. (2009), but at the same time there are 
also several differences. For instance, with respect to the refer- 
ence PPM entropy profile, pre-shock entropy production and a 
lower entropy level behind the shock is common to both of the 
tests. However, in Wetzstein et al. (2009), pre-shock heating for 
the time-dependent AV models occurs at radii larger than for the 



12 



R. Valdarnini: The impact of numerical viscosity in SPH simulations of galaxy clusters 



0.001 



- 0.0001 
II 

io- 5 



o 



10- 



io- 



io- 



0.001 



~ 0.0001 

II 

<: 10- 5 



io- B 
io- 7 

IO" 8 




-\ — I — h 




10 

K kp /(2 tt) 



10 

K kp /(2 tt) 



Fig. 5. Compressive components of the density-weighted velocity power spectra (20) are shown at z = as a function of the dimensionless 
wavenumber k = k r L sp /2n, where k r = \k\. The spectra are extracted from the high resolution (HR) runs of the four dynamically perturbed test 
clusters using a cube of size L sp = r 2 oo with = 128 3 grid points and are shown up to the wavenumber k = N g /2. In each panel, different lines 
refer to runs with different AV viscosity parameters, the line style and color coding is the same as in the previous figures. The solid black line 
indicates the Kolgomorov scaling, while open diamonds correspond to the volume-weighted velocity power spectrum of the standard AV runs. 



fixed AV run and the differences in the entropy profiles (as well 
as in the density) for different AV runs are here reduced more. 

It is difficult to ascertain the origin of the different behav- 
ior of the two codes in the proximity of the shock front, given 
also the fact that the differences in the corresponding profiles are 
not dramatic. However, the SPH formulation implemented in the 
two codes differs in several aspects, so that even results produced 
for the same test problem and using the same initial conditions 
might be different. The thermal evolution of the gas is followed 
here according to an entropy-conserving scheme, while in the 
Wetzstein et al. (2009) code it is the specific internal energy 
which is integrated. Moreover, unlike Wetzstein et al. (2009), 
the equation of motion (3) incorporates the terms due to the 
presence of smoothing length gradients. The reader is referred 
to Springel & Hernquist (2002) for a thorough discussion of the 
relative differences between simulation results produced using 
the two SPH formulations. 

To summarize, for the test problem considered the results 
presented in this section agree with previous findings and vali- 
date the code, as well as showing its capability to properly model 
shocks which develop during the collapse of self-gravitating ob- 
jects. In accordance with the results of Sect. 5.1.2, profiles ex- 
tracted from simulations with a time-dependent AV scheme ex- 



hibit a behavior which is very similar to the corresponding ones 
for the fixed AV model, thus showing shock resolution capabil- 
ities analogous to those of the standard AV scheme for these 
models. Finally, the peaks of the viscosity parameters at the 
shock front appear almost independent of the chosen value of 
the AV decay parameter Id, thereby supporting the proposed 
parametrization (30). 



6. Cluster simulations 

In this section the impact of numerical viscosity on the ICM ve- 
locity field of simulated clusters is studied using the statistical 
tools presented in Sect. 4. The purpose of the analysis is also to 
obtain useful indications about the numerical constraints which 
are necessary in order to adequately describe turbulence in sim- 
ulations of the ICM using SPH codes. 

6. 1 . Velocity power spectra 

Fourier spectra extracted from hydrodynamical simulations are 
expected to exhibit a dependence on both the simulation reso- 
lution and the hydrodynamical method used in the simulations. 
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Fig. 6. The same as in Fig. 5, but for the shearing components of the velocity power spectra. 



Usually, resolution issues are debated in a dedicated section with 
the purpose of assessing the stability of the simulation results. 
Here, however, the interpretation of the results is strictly related 
to the resolution employed in the simulations, so that the two 
arguments will be discussed together. 

In order to measure the velocity Fourier spectra of the sim- 
ulated clusters, a cube of side length L sp is placed at the cluster 
center, the latter is defined as the location where the gas den- 
sity reaches its maximum value. As described in Sect. 4, Fourier 
transforms of the density weighted velocity field are evaluated 
by first estimating interpolated quantities at the cube grid points 
and then performing a 3D FFT of the sampled values. The choice 
of the cube side length L sp and of the number of grid points iVj is 
however limited by several arguments which strongly constrain 
the possible choices. 

At variance with studies of supersonic turbulence using hy- 
drodynamical simulations (Kitsionas et al. 2009, herefater K09; 
Price & Federrath 2010) here the gas distribution is driven by 
gravity and because of the Lagrangian nature of the SPH code, 
the bulk of the particle distribution is located in the central clus- 
ter regions. For the cluster simulations presented here, at the 
present epoch about half of the cluster mass is contained within a 
radius of ~ 7"2oo/3. Therefore, in order to minimize this aspect of 
the resolution effects, the size of the cube should be kept as small 
as possible, but in this case most of the large-scale modes will 
be missed in the spectral analysis, in particular, those modes cor- 



responding to the merging and accretion processes of the clus- 
ter substructure. As a compromise between these two opposite 
needs, the side-length of the cube is set to L sp - r2oo (the scaling 
L sp 00 ^"200 is chosen in order to consistently compare velocity 
spectra extracted from different clusters). 

The choice of the number of grid points 2V? in the cube or , 
equivalently, of the grid spacing L sp /N g , depends on the effec- 
tive SPH resolution of the simulations. In SPH the smallest spa- 
tially resolvable scale is set by the values of the gas smoothing 
lengths hi. As already outlined, because of the Lagrangian nature 
of SPH simulations, the gas particle number density increases in 
high-density regions and this in turn implies a subsequent de- 
crease in the gas smoothing lengths hi. If in these regions the 
SPH spatial resolution becomes higher than the grid resolution 
of the cube, SPH kernel interpolation at the cube grid points will 
appear in the fc-space as a extra power at the highest wave num- 
ber permitted by the grid. This effect has been noticed by K09 
and can be avoided provided that the cube grid spacing is chosen 
suitably small. For the HR cluster simulations, values of the gas 
smoothing lengths h, in the cluster cores range from ~ 20 kpc 
for the most massive clusters down to hi ~ 5 kpc for the least 
massive ones. Similar values for the grid spacing are obtained 
setting N g = 128 grid points along each of the spatial axis of the 
cube. However, the number of cube grid points cannot be made 
arbitrarily high because of the need for avoiding undersampling 
effects in the estimate of SPH variables at the grid points. As a 
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rule of thumb, the number N g should not exceed ~ 2N p , with 
N p being the number of SPH particles. 

Finally, when calculating the velocity power spectra non- 
periodicity effects should be taken into account by adopting a 
zero-padding technique (Vazza et al. 2009). However, applica- 
tion of this method to the clusters presented here requires some 
care since the choice of the cube side length L sp is dictated by the 
previous constraints so that a cube of volume L^ sp occupies at the 
cluster center ~ 1 /4 of the volume of a sphere of radius r^oo with 
its origin centered as the cube. Hence, a power estimate based 
on the non-periodicity assumption would likely underestimate 
the velocity power spectrum at wavenumbers ~ nlL sp , in par- 
ticular for those clusters which are undergoing a major merging 
event. The procedure adopted here should then provide a better 
indication of the spectral behavior of the ICM velocity field at 
spatial scales of the order of the cube size. 

Having chosen the cube parameters, velocity power spectra 
have been evaluated at the present epoch for the ensemble of 
simulated clusters constructed by collecting simulations of the 
the cluster HR baseline sample performed with different AV pre- 
scriptions. The resulting spectra are shown in Figures 5 to 8 and 
provide a quantitative comparison between statistical properties 
of the longitudinal and solenoidal velocity field spectral compo- 
nents for clusters simulated with different AV settings and with 
different dynamical states. Each panel in the figures corresponds 
to an individual cluster and the curves shown in the panels are 
the velocity power spectra extracted from simulations of the con- 
sidered cluster realized using different AV parameters. 

A first conclusion to be drawn from the spectra shown in 
the figures is that, for a given cluster and power spectrum com- 
ponent, at large scales the velocity power spectra do not show 
systematic differences between different AV runs and therefore 
the effects of numerical viscosity can be considered negligible 
on these scales. Moreover, dissipative effects in the kinetic en- 
ergy are relatively less important for those runs in which the pa- 
rameter settings of the time-dependent AV scheme correspond 
to the shortest decay time scales for the viscosity parameter a(t). 
In fact, at high wavenumbers the power spectra of the AV 5 runs 
are higher than those of the standard AV runs (AV ) by a factor 
of ~ 2 and in certain cases even by a factor of ~ 10 ( see, for 
example, the longitudinal power spectrum of cluster 110 of the 
relaxed subsample). These behaviors are shared by the velocity 
power spectrum components of all of the simulated clusters and 
because of the sample size can be considered as systematic. 

The use of Fourier spectra allows one to study in a quanti- 
tative way the scale dependency of dissipative effects due to nu- 
merical viscosity. All of the spectra exhibit a maximum at k < 2, 
at spatial scales of ~ f20u/2. This is a key signature in the power 
spectra of the ICM velocity field, in which merging and accre- 
tion processes driven by gravity at cluster scales are the primary 
sources of energy injection into the ICM. Similar findings have 
been obtained by Vazza et al. (2009), who measured the veloc- 
ity power spectrum of simulated clusters in cosmological simu- 
lations using the AMR Eulerian code ENZO. 

Spectra of runs with different AV settings begin to deviate at 
wavenumbers k > 6 = kj. This indicates that for the present sim- 
ulation resolution (HR) the effects of numerical viscosity on the 
formation of eddies due to Kelvin-Helmholtz instabilities gen- 
erated by substructure motion (Takizawa 2005; Subramanian, 
Shukurov & Haugen 2006) become significant at spatial scales 
Idiss = x/kd ~ ?"2oo/10 ~ 100 - 200 kpc, in accordance with the 
behavior of the velocity structure functions (see later). The gen- 
eration of these instabilities leads to the development of turbu- 
lent motions which, as expected, are stronger in those runs with 



the lowest numerical viscosity. There are no large differences 
seen between the shape of the spectra extracted from clusters 
in different dynamical states. However, the power spectra of re- 
laxed clusters, in comparison with the spectra of the dynamically 
perturbed clusters, show a certain tendency to flattening at large 
scales and smaller amplitudes at the smallest cube wavenumber. 

For comparative purposes, each panel shows for the con- 
sidered case the volume-weighted velocity power spectrum 
of the standard AV run. The behavior of the other volume- 
weighted spectra reflects the similarities seen in the correspond- 
ing density-weighted spectra of simulations with different AV 
settings and for the sake of clarity these spectra are not shown in 
the panels. 

At large scales, volume-weighted and density-weighted 
spectra do not show significant differences, with few expec- 
tions such as the longitudinal power spectra of cluster 105 of 
the perturbed subsample. Variations in the behavior of density- 
weighted spectra and volume-weighted spectra are expected to 
appear at smaller scales, when gas density gradients become sig- 
nificant. Note, however, that there are no systematic differences 
between density and volume-weighted spectra as a function of 
wavenumber versus the cluster dynamical state and/or power 
spectrum decomposition. 

At high wavenumbers, in certain cases density-weighted 
spectra can have values below those of volume-weighted spectra. 
This is interpreted as a resolution effect, which is stronger as runs 
with lower resolution (MR or LR) are considered. As already 
outlined, the gas distribution is driven by gravity and because of 
the Lagrangian nature of the SPH code, for wavelength modes 
approaching the grid spacing the contribution of low density re- 
gions which are undersampled becomes dominant. This effect is 
not present in all of the clusters because it is strongly dependent 
on the degree of gas inhomogeneity which in turn depends on 
the individual cluster dynamical history. 

Although the aim of the paper is to assess the effectiveness 
of the time-dependent AV scheme in reducing the numerical 
viscosity in SPH cluster simulations, nevertheless the velocity 
power spectra presented here can be used to obtain useful in- 
dications about the description of turbulence in the ICM when 
SPH simulations are used. 

Below the driving scales k ~ 1-2, the longitudinal com- 
ponent of the energy spectra exhibit an approximate power- 
law behavior close to the Kolgomorov scaling that extends 
down to dissipative scales. These features are not shared by the 
solenoidal component of the spectra. To study the behavior of 
the shearing part of the spectra, it is useful to introduce the ratio 
^(k) = E c (k)/(E c (k) + E s (k)) which measures as a function of 
the wavenumber the fractional contribution of the longitudinal 
velocity power spectrum component. This ratio lies in the range 
^(k) ~ 0.3 - 0.5 at large scales and rises above k ~ 5 = k r s es , 
the latter is defined as approximately the wavenumber at which 
occurs the bending in ^(k). While the ratio at large scales can 
be understood in terms of fully developed turbulence, the rise 
in ^(k) at small wavenumbers is a consequence of the steep de- 
crease in the shearing part of the spectra in this range of scales, 
which is interpreted here as a resolution effect. This behavior 
shows that a fundamental issue to be addressed is the resolu- 
tion dependence of the velocity power spectra extracted from 
the simulations. 

K09 and Price & Federrath (2010) studied supersonic 
isothermal turbulence using both SPH and Eulerian based codes. 
Both groups reach the conclusion that differences in the simula- 
tion results are mainly due to the resolution used, with dissipa- 
tive effects specific to the codes confined at high wavenumbers. 
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Fig. 7. The same as in Fig. 5, but the spectra here are computed for the quiescent test clusters. 



In fact, equivalent results in the measured velocity power spectra 
are obtained using an equivalent number of resolution elements 
N for each direction in space in simulations with different codes. 
For example, SPH simulations with A^ 3 = 512 3 particles would 
give comparable results to those obtained using the same number 
of grid cells in Eulerian simulations. 

Exploiting this correspondence one can then interpret the 
dependence of the velocity power spectra on resolution effects 
using the results of Federrath et al. (2010), who studied inter- 
stellar turbulence using a FLASH code with up to A^ 3 = 1024 3 
cell elements. From the simulations the authors (see their sect. 
6) measure in the extracted Fourier spectra E(k) a rise in *?(&) 
starting at k ~ 40 ~ N/26 which they explain as a resolution 
effect. Basically, this happens because the number of computa- 
tional elements needed to accurately resolve rotational modes 
is larger than that necessary to describe longitudinal modes, for 
which only one degree of freedom is present. In fact, longitu- 
dinal modes appear to be well described down to k ~ N/10, 
suggesting a difference in the resolution scale by about a factor 
of ~ three. 

Application of these results to the present simulations sug- 
gests that the solenoidal part of the spectra is barely resolved in 
the HR runs, with a resolution scale k r s es ~ N/26 ~ 64/26 ~ 3, 
whilst the longitudinal component is resolved down to k r c es ~ 
N/10 ~ 6. To corroborate these findings, Fig. 9 shows the 
solenoidal component of the density-weighted velocity power 



spectra versus the wavenumber for simulations of different res- 
olution for the subsample of relaxed clusters. All of the spectra 
show a well-defined tendency towards shallower slopes as the 
simulation resolution is increased, with a reduction of the dissi- 
pative range towards higher wavenumbers. This spectral behav- 
ior is also partly explained by the larger number of matter clumps 
which are now resolved at small scales as the mass resolution is 
increased, thus adding power to the spectra at high wavenum- 
bers. At large scales a slight decrease in the amplitude of the 
spectra can occur because of resolution effects, as the resolution 
is increased clumps of matter accreting at r ~ rajo survive more 
easily therefore reducing the amount of turbulence injected into 
the medium via ram pressure stripping mechanisms. This effect 
can be seen in Fig. 9 for clusters Oil and 019, for which the large 
scale part of the spectra decreases as the resolution is increased. 
Moreover, note that in certain cases there is a plateau in the spec- 
tra at high wavenumbers which is removed as the resolution is 
increased. 

Fig. 10 shows the ratio ^(k), for the same test runs. Although 
there are large cluster-to-cluster variations, nevertheless cluster 
110 offers a particularly neat example of the shift of k' s cs towards 
higher wavenumbers as the resolution is increased. This depen- 
dence on the numerical resolution also confirms that the pre- 
dominance of the compressive components in the velocity power 
spectra in not due to the presence of shocks at small scales, with 
the primary injection mechanisms occurring at cluster scales. 
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However, the estimates of k' s es inferred from the application of 
other simulation results might be too pessimistic. From Fig. 10 
it can be seen that for clusters 013 and 110, ^(k) is stable down 
to k ~ 6 - 7. These values are a factor of ~two higher than those 
estimated from previous scalings. 

Similar underestimates are found for the dependence of the 
dissipative scale set by numerical viscosity versus the numerical 
resolution. According to K09, length-scales defined by the con- 
dition k > N/32 are affected by dissipative mechanisms specific 
to the code. This implies that dissipative effects here should start 
acting from k ~ 2, whereas the spectral behavior indicates that 
at least down to k ~ 6 = kd there are no significant differences 
between spectra extracted from runs with different AV settings. 

An important issue which needs to be addressed is the de- 
pendence on the numerical resolution of the dissipative scales 
identified from the spectral behavior of the velocity power spec- 
tra. For doing this the use of the longitudinal power spectrum 
components is more suitable because for these spectra, numeri- 
cal resolution effects are less severe than for the solenoidal com- 
ponents. In fact, a comparison as in Fig. 9 of the behavior of the 
longitudinal part of the density-weighted velocity power spec- 
tra for simulations of different resolutions, shows that the dis- 
sipative wavenumber kd scales approximately as N, suggesting 
that in order to resolve eddies down to length-scales of the or- 
der of galaxy-sized objects (~ 50 kpc) simulations would be 
required with at least A^ 3 = 256 3 gas particles. However, the 



spatial range of l t a ss derived from the spectra of Figures 5-8 is 
in rough agreement with that expected using analytical models 
(Subramanian, Shukurov & Haugen 2006) to estimate the co- 
herence velocity field length-scale associated with minor merg- 
ers in clusters. Moreover, Maier et al. (2009) have recently pre- 
sented an Eulerian-based AMR code with subgrid modeling for 
turbulence at unresolved scales. The authors studied the evolu- 
tion of turbulence in cluster simulations and argue that most of 
the turbulent energy is injected into the medium at length scales 
of ~ 125 - 250/T kpc (cf. their Fig. 6). This consistency be- 
tween different independent estimates of the characteristic turbu- 
lent length scales therefore suggests that for the HR runs, spec- 
tral features due to dissipative effects might be of physical origin 
rather than due to numerical resolution effects. 

To summarize, application of previous findings (Kitsionas et 
al. 2009; Federrath et al. 2010; Price & Federrath 2010) to the 
simulations presented here produces several discrepancies in the 
resolution dependence of characteristic scales, which appear less 
severe here than those found by these authors. A clarification of 
these issues clearly has to await simulations with much higher 
resolution; nevertheless it is suggested that physical differences 
in the problems considered might play a role. For instance, grav- 
ity is here the driving force which dominates the dynamics of the 
collapsed object, thereby creating, owing to the adaptative nature 
of the SPH codes in space, a 'scale dependence' in resolution as 
one moves outwards from the cluster center. As a consequence, 
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Fig. 9. Shearing components of the density-weighted velocity power spectra are shown at z = for the runs of the relaxed test clusters as a 
function of the wavenumber. As in the previous figures, different line styles are for different AV runs. In each panel runs with different resolution 
are displayed with different line thickness. For the sake of clarity, not all of the AV runs are plotted and runs with different AV parameters have 
been equispaced from top to bottom by a factor of 10. To better illustrate the wavenumber dependence, spectra have been computed here using 
Nji = 25 6 3 grid points. 



the effective resolution of the simulations performed here might 
be higher than that estimated using simple scaling arguments. 

However, these results indicate that the present simulation 
resolution is clearly inadequate to fully describe the spectral be- 
havior of the ICM velocity field, therefore raising the question of 
what is the minimum number of particles needed in order to ad- 
equately resolve turbulence in SPH simulations. This is strictly 
related to the minimum spatial scale which must be resolved in 
the velocity Fourier spectra. Requiring for rotational motion at 
least one order of magnitude of spectral resolution, implies that 
N 3 = 512 3 particles are needed for a correct description of tur- 
bulence in SPH simulations of the ICM. Taking into account 
the previous caveats, the same results could be obtained using 
N 3 = 256 3 particles. 

In addition to spectral analysis, the scale dependence of dis- 
sipative effects can be studied using the velocity structure func- 
tions which provide information in physical space about the ve- 
locity field self-correlation. In particular, the second order struc- 
ture function Szff) is complementary to Figures 5-8 for obtain- 
ing informations about the effects of numerical viscosity on the 
velocity field statistics. To evaluate the structure functions is a 
computationally demanding task since the number of pairs scales 
approximately as oc N? as . The choice of the subsample size N s is 



a compromise between the needs of adequately sample the struc- 
ture functions and of reducing the computational cost. The par- 
allel and transverse second order structure functions have been 
calculated here setting N s = N gas /8, corresponding approxi- 
mately to ~ 5 • 10 9 pairs for the HR runs. This choice for N s 
is found to give convergent results when tested against higher 
values. Note that a uniform random number extraction of N s 
particles from the simulation particles does not correspond to 
an approximately uniform sampling in space, as the bulk of the 
gas particles (~ 1 /2) are concentrated within distances i f"2oo/3 
from the cluster center. 

The parallel and transverse second order velocity structure 
functions are shown in Fig. 1 1 for the HR runs using volume- 
weighted velocities. All of the functions exhibit a growing be- 
havior with increasing pair separation f = r/r2oo, thereby show- 
ing how in the ICM turbulence is substained by motions on clus- 
ter scales. For two of the clusters (1 and 5), the radial depen- 
dence is approximately a power-law, a result which can be taken 
as indicative of a scaling range, although this issue needs to be 
clarified with simulations of much higher resolution. 

Assuming as reference profile that of the structure functions 
produced by the standard AV runs, structure function profiles 
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Fig. 10. The ratio of the longitudinal to total velocity power spectra are shown at z = for runs of the relaxed test clusters as a function of the 
wavenumber. The line coding is the same as in the previous figure; runs with different AV parameters have been equispaced from top to bottom by 
a factor 0.2. 



extracted from the runs with a time-dependent AV scheme at 
small scales deviate significantly from the reference profile. The 
deviation increases as f decreases and is wider as the decaying 
viscosity parameter lj approaches unity. This behavior clearly 
shows the effectiveness of the new numerical viscosity scheme 
in reducing the viscous damping of velocity away from shocks, 
and therefore of being the code better able to follow the devel- 
opment of turbulent motion. Profiles of the structure functions 
extracted from simulations with lower resolution are noisier than 
the profiles presented here but preserve the basic features shown 
in Fig. 11, therefore supporting the validity of these arguments. 



Structure functions derived using a density-weighted veloc- 
ity are shown in Fig. 12. They exhibit the same behavior as for 
the volume-weighted functions but with a much shallower ra- 
dial dependence. This is due to the adopted weighting scheme, 
for which the contribution of those pairs located in high den- 
sity regions is enhanced with respect to the volume-weighted 
functions. This preferentially happens at small scales because in 
SPH codes most of the particles are located in high density re- 
gions. Finally, a comparison with the spectra shown in Figures 
5-8, shows an approximate correspondence between the length 
scales (~ r2oo/10) at which dissipative effects become signifi- 
cant. 



6.2. Average radial profiles 

Measurements of the average radial profiles of thermodynamic 
variables, such as density or temperature, extracted from spa- 
tially resolved X-ray observations allow one to probe ICM prop- 
erties and test the predictions of cosmological cluster simula- 
tions. This is achieved by comparing the observed ICM profiles 
with those extracted from mock X-ray data. It is therefore impor- 
tant to examine in SPH cluster simulations the impact of numer- 
ical viscosity on the radial profiles of the variables considered. 

For doing this, the final radial profiles of the gas entropy 
S(r) = k.BT(r)/[im p p^ 3 are shown in Fig. 13 for the simulated 
clusters of the relaxed subsample. As in Figures 5-8, each panel 
corresponds to an individual cluster and within each panel dif- 
ferent curves are for runs with different AV settings. In order 
to compare the entropy profiles of different clusters, the entropy 
has been normalized to 
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where fi, = Q.^/Q. m is the global baryon fraction. 

Discussion of the dependence of the final profiles on the AV 
scheme used in the simulations is restricted here to the entropy 
only, since both density and temperature profiles exhibit quite 
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Fig. 11. Transverse (5]|) and longitudinal (S±) second-order normalized velocity structure functions at z = for HR runs of the perturbed test 
clusters as a function of r/r 2 oo- 



similar behavior. A feature common to all of the entropy profiles 
is the relatively small scatter between those profiles extracted 
from runs for a given cluster with different strengths of the nu- 
merical viscosity. The variations in the entropy profiles due to 
the AV implementation appear to be confined to the inner cluster 
regions (~ 10% of r2oo)- As already outlined in the Introduction, 
a fundamental issue to be investigated is the amount of entropy 
present in cluster cores (Agertz et al. 2007; Wadsley, Veeravalli 
& Couchman 2008; Mitchell et al. 2009). Viscous damping of 
random gas motion is expected to be reduced in low-viscosity 
runs, thereby increasing fluid mixing and possibly the accretion 
at the cluster center of gas with higher entropy, with a subsequent 
rise of the entropy floor in the core. 

The profiles of Fig. 13 illustrate that the amount of gas mix- 
ing due to numerical viscosity is not significant, although there 
is a certain tendency for the core entropy to be influenced by the 
AV strength of the simulation. For instance, cluster 1 1 exhibits a 
large difference of ~ 5 between the core entropy S (r c = 0.01r2oo) 
of run AVo and that of run AV4. For two other clusters (12 and 
1 10), the entropy profiles of the low-viscosity runs have a mod- 
erate increase with respect to that of the standard AV run as 
the radial coordinate decreases, with central differences being 
smaller than a factor of ~ two. Finally, there is one cluster (19) 
for which S(r c ) of run AV4 is below that of run AVo. This be- 
havior is shared by the entropy profiles of the perturbed cluster 



subsample and it indicates that the amount of entropy mixing in 
cluster cores induced by numerical viscosity effects is modest or 
quite negligible. 

The stability of the results is assessed in Fig. 14 by contrast- 
ing for the same test clusters the entropy profiles of simulations 
with different resolution. In each panel the profiles of only three 
AV runs are plotted (AV2, AV4 and AV5) to avoid overcrowding. 
From Fig. 14 it can be seen that there is little variation between 
the entropy profiles of runs performed with the same AV set- 
tings but with different resolution. In the AV5 test case, for sev- 
eral clusters the profiles of the HR runs show some deviations at 
small distances from those with lower resolution. Interpretation 
of this difference as a systematic effect due to inadequate resolu- 
tion requires some care however, since it is absent in the profiles 
of the other AV runs. Moreover, even for the test case AV5, in 
one cluster (110) the core entropy S(r c ) does not vary with the 
simulation resolution. 

To summarize, the reliability of the entropy profiles of sim- 
ulated clusters with respect to variations in the simulation reso- 
lution is supported by the profiles of Fig. 14, for which one sees 
the lack of any systematic dependence of S (r e ) with varying nu- 
merical resolution. This indicates that the profiles of Fig. 13 can 
be reliably used to draw the conclusion that in SPH simulations 
the amount of entropy present in cluster cores due to numerical 
viscosity effects can be considered negligible. The smallness of 



20 



R. Valdarnini: The impact of numerical viscosity in SPH simulations of galaxy clusters 



II 



in 



=C\2 

m 



0.1 



0.01 




0.001 

1 



0.1 



m 



in 



0.01 



0.001 



~l — I — III 1 1 — I — I — III 



001 



PHR 



I I I I L I 



H 1 I I I I I I 




016 



_i i i 1 1 1 1 

0.1 



_l I I I I I I I 



~l — I — I — I I I 



~l 1 1 — I — I I I I 



005 



105 



_i i i i 1 1 1 1 

0.1 



_J I I I I I I I 



r/ r 



200 
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AV effects in determining the core entropy of simulated clusters 
therefore indicates that the discrepancy, which originates at the 
cluster centers between SPH and grid-based Eulerian codes in 
the produced entropy profiles of adiabatic simulations, is inher- 
ent to the method itself rather than being a problem of the AV 
scheme implemented in the SPH codes. In particular, the lack of 
significant AV effects indirectly confirms that the primary source 
of the entropy discrepancy between the two hydrodynamical ap- 
proaches resides in the treatment of flow discontinuities in SPH 
simulations (Agertz et al. 2007). 

To overcome these difficulties Price (2008) and Wadsley, 
Veeravalli & Couchman (2008) proposed to add an explicit dif- 
fusion term in the SPH thermal energy equation, with the pur- 
pose of smoothing the discontinuities present in the thermal en- 
ergy at the fluid interfaces. The introduction of this diffusion 
term leads to a certain amount of gas mixing and in turn to an 
establishment of an entropy core in simulated clusters (Wadsley, 
Veeravalli & Couchman 2008), which appears to qualitatively 
agree with that produced by Eulerian codes. However, genera- 
tion of vorticity associated with fluid instabilities is severely af- 
fected by viscous damping due to numerical viscosity. Using a 
suite of 2-D numerical tests Price (2008) showed that the treat- 
ment of Kelvin-Helmholtz instabilities can be greatly improved 
in SPH codes provided that the hydrodynamic equations are re- 
formulated so as to include thermal diffusion as well as variable 
AV coefficients. It would then be interesting to reinvestigate for 



the same test sample presented here the growth of central en- 
tropies using SPH simulations in which the hydrodynamic equa- 
tions have been modified by adding an artificial thermal energy 
dissipation term. 



The radial profiles a(r) of the viscosity parameters are shown 
in Fig. 15 for the relaxed subsample. For each cluster, only the 
profiles of runs AV2 and AV5 are shown, the radial behavior of 
the profiles for the other AV runs being intermediate between 
the two. Within each panel, the profiles of the AV test cases are 
accompanied by the viscosity profiles extracted from the corre- 
sponding low and medium resolution runs. This is in order to 
estimate the stability of the profiles with respect to varying the 
numerical resolution of the simulations. From the figure it can 
be seen that for a given AV setting there is little scatter between 
the profiles extracted from runs with different resolutions, with 
some differences near to the shock front which are not signif- 
icant however. The profiles of the perturbed cluster subsample 
are similar to those shown here, but with a wider scatter because 
of the different dynamical history of these clusters. The radial 
behavior of the profiles is in agreement with what was expected: 
note that for two clusters (19 and 1 10), the peak values apeak of 
a(r) at the shock front for the run AV5 are significant smaller 
than those for run AV2. This is not the case for the peak values 
of run AV4 (not shown here), which are similar to those of run 
AV2, indicating therefore that for weak shocks (say a pea k i 0.3) 
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Fig. 13. Final radial entropy profiles for the HR runs of the relaxed test clusters. Gas entropy is defined as S(r) = k B T(r)/ fim p pg and is plotted 
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the peak value of the viscosity parameter is influenced by the 
floor value a ml „ of the gas particles ahead of the shock. 

One of the most important aspects in which viscous damp- 
ing of random gas motion is expected to influence ICM prop- 
erties, is in the ICM energy budget. In order to investigate the 
impact of numerical viscosity on the energy content of the ICM, 
it is necessary however to properly disentangle from the velocity 
field u(x) the contribution of bulk motions from the small-scale 
chaotic parts. Following Dolag et al. (2005), a turbulent velocity 
field u(x) can then be defined by subtracting from u(x) a local 
mean velocity < u(x) >. According to Eq. 25, the latter is de- 
fined by using a TSC window function to interpolate gas density 
and velocities at the grid points of a cubic mesh with grid spac- 
ing H. The mean velocity < u(x) > is defined by averaging over 
the gas particles which belong to the grid cell; subtraction of the 
filtered velocity from the velocity field u(x) is expected to re- 
move the contribution of large scale motions and leave a velocity 
u{x) with spectral content dominated by those wavenumbers for 
which kH » 1. The choice of the filtering scale H is somewhat 
arbitrary: Dolag et al. (2005) found that a mesh spacing of the 
order of ~ 30 kpc yields consistent results; in agreement with 
the choices of Sect. 6.1 the scaling H oc r2oo is adopted here, 
allowing the velocity field u(x) of clusters of different masses to 
be compared in a self-similar way. The choice of H is a com- 



promise between the needs of removing the contribution of bulk 
flow motions to u(x) as much as possible and the constraints im- 
posed by the numerical resolution of the simulation for which 
values of H which are too small leads to an undersampling of 
local estimates. These effects are expected to be significantly re- 
duced by requiring at least N c ~ 10 2 gas particles in the cube 
cells. From this constraint, the allowed range for the grid spacing 
H can be estimated using the average gas density radial profiles. 
For the clusters presented here, p/p c ~ 5 at the cluster radius 
r2oo and for the HR runs this translates into a range of values 
for H between ~ 740 kpc for the most massive clusters down to 
~ 350 kpc for the least massive ones. 

Such filtering scales are too large to effectively remove bulk 
motion components from the velocity field u(x). As a compro- 
mise, the contribution of turbulent energy can be investigated by 
restricting the study to the inner parts of the cluster. At the radius 
r ~ ?"2oo/2 the cluster gas densities lie around p/p c ~ 30 and the 
range for H is now between ~ 400 kpc and ~ 180 kpc, which 
corresponds to H ~ r2oo/5, when rescaled to the cluster radius 
r 200- Although this choice does not completely remove from ii(x) 
the laminar flow patterns present in the ICM, nevertheless it is 
supposed to filter out most of the contribution to u(x). In a recent 
paper Vazza et al. (2009) studied the development of turbulence 
in the ICM using the adaptative grid AMR cosmological code 
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Fig. 14. As in Fig. 9, but for the entropy profiles. 



ENZO and argued that for their simulated cluster the ICM ve- 
locity field at length scales i 300 kpc is dominated by laminar 
flows. Moreover, the turbulent energy density radial profiles con- 
structed from a filtered field using a length scale of H/2 do not 
show appreciable differences from those derived here using H as 
a filtering scale. 



In order to consistently compare our simulation results with 
those of Vazza et al. (2009), the radial density profiles are con- 
structed using the same definitions. The thermal energy den- 
sity is defined as E,i, = 3kBT{r)p g (r)/2pm p , where T(r) is 
the mass-weighted gas temperature, kg is the Boltzmann con- 
stant, jj = 0.6 is the molecular weight and m p the proton mass; 
Ekin = p g u 2 /2 and E tur b = p g u 2 /2 are, respectively, the kinetic 
and turbulent energy density. The total energy density is defined 
as E,„i = E,h+Ekm- All of the quantities with a radial dependence 
are estimated by averaging over a collection of values; these are 
calculated according to the SPH prescription (19) at a set of grid 
points uniformly spaced in angular coordinates and lying at the 
surface of a spherical shell, located at distance r from the clus- 
ter center. The energy density profiles constructed in this way 
are shown in Fig. 16 for the clusters of the perturbed subsam- 
ple. The profiles have been rescaled in units of p c cr 2 and for the 
sake of clarity for each cluster only the profiles extracted from 
the runs of two AV test cases (AVo and AV4) are shown. The 
corresponding ratios E tur b/E m are shown in Fig. 17. 



A first conclusion to be drawn from the radial behavior of 
the profiles is the significant impact on the ratio E tur b/E to , of 
the numerical viscosity scheme adopted in the simulation. As a 
consequence of the suppression of viscous damping of random 
gas motion, for the runs AV4 the ratio E tur b/E to , at any radial 
distance is higher than that of the corresponding runs AVo by a 
factor ranging between ~ 30% and ~ 100%. The ratio is nearly 
constant for the cluster 005 ( ~ 8%) and 016 ( ~ 10% ) and 
increases from ~ 2% to ~ 5% for cluster 105. Finally, for cluster 
001 it decreases (~ 10% at r — r2oo/100) as r increases (~ 5% at 
r = r2oo/2). A similar behavior is present for the profiles of the 
relaxed cluster subsample. 

These findings indicate the lack of any systematic depen- 
dence of the ratio E tur b/E tot on the cluster radial distance and 
illustrate the sensitivity of the turbulent energy density profile 
on the dynamical history of the individual cluster. In order to 
extract statistically meaningful results about the contribution of 
random gas motion to the ICM energy budget, it is then neces- 
sary to use a very large simulated sample ( say £ 100 clusters, 
see, e.g., Piffaretti & Valdarnini 2008). 

It is difficult to assess the consistency of the profiles pre- 
sented here with those obtained by Vazza et al. (2009), since 
the results indicate large cluster-to- cluster variations among the 
individual turbulent energy density profiles and in Vazza et al. 
(2009) the profile is shown for only one cluster. However, the 
authors argue that, for their cluster, the turbulent energy budget 
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Fig. 15. As in Fig. 9, final profiles of the viscosity parameter a are shown for the runs of the relaxed test clusters with different AV parameters and 
resolution. In each cluster panel the profiles have been plotted without any shifting between the various cases. 



is smaller by a factor ~ 5 - 6 than previous SPH results. The 
ratio E, ur b/E to , lies here in the range between a few percent and 
~ 5 — 10%, and similar values for this ratio are quoted by Vazza 
et al. (2009). Moreover, a visual inspection shows a strong sim- 
ilarity between the profiles of cluster 105 and those derived by 
Vazza et al. (2009) for their simulated cluster ( Fig. 6, bottom 
right). As already stressed in Sect. 6.1, a proper study of tur- 
bulence in the ICM using an SPH code awaits simulations of 
much higher resolution than those presented here, nevertheless 
these results are encouraging and support the view that SPH sim- 
ulations can be used to investigate this topic, provided that the 
numerical viscosity effects are properly treated and the medium 
is adequately resolved. 

An important issue which needs to be investigated is the im- 
pact of the numerical viscosity scheme on the cluster mass bias 
estimated from simulations. Recently, the relative importance 
of different mass bias components has been studied in detail 
(Piffaretti & Valdarnini 2008) using a large set of N-body/SPH 
simulations of galaxy clusters. Apart from the bias originating 
from the use of observational quantities, such as spectroscopic 
temperatures, and those induced by assuming specific models 
for the measured profiles, it was found that mass underestimates 
are dominated by the hydrostatic equilibrium assumption and by 
the absence in the Jeans equation for the mass of pressure sup- 
port terms due to random gas motions. In particular, for relaxed 



clusters the mass bias due to these terms increases with the dis- 
tance r from the cluster centre and can be as high as ~ 5 - 10% 
for r ranging between r25oo and rjjoo. Given that the importance 
of galaxy clusters as cosmological probes relies on the use of 
accurate mass estimators, it is therefore important to estimate in 
the Jeans equation the dependency of these gas motion terms on 
numerical viscosity. 

For doing this, the relative importance of the velocity pres- 
sure terms is quantified by constructing the radial profile of 



Pa 



Pa + Pa 



cr 2 (r)/3+kT(r)/iJm p 



(33) 



where <x(r) is the rms gas velocity dispersion. The latter 
is evaluated at the same radial shells used to evaluate the en- 
ergy density profiles. In analogy with Eq. 33, a relative turbulent 
pressure term can be defined, P tU rb/Ptou f° r which Pfurb/Pm = 
Eturb/Etot- The quantities Pa-/(P<r+Pth) can then be directly com- 
pared with the turbulent-to-total energy density ratios and these 
are shown in Fig. 17 as thin lines. For the four test clusters, the 
ratio increase with the radial distance r and Pa- 1 {Pa + Pth) ~ 
10 - 20% at r = rsno ~ 0.6r2oo> while at the same distance 
Pa/(Pa + Pth) ~ 10% for the relaxed subsample. 

This range of values for the contribution of random gas mo- 
tions to the total pressure is in accordance with the estimates 
extracted from previous simulations (Rasia et al. 2006; Nagai, 
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Fig. 16. Final energy density radial profiles constructed from the HR runs of the dynamically perturbed test clusters. Each panel shows: the profile 



of the thermal energy density E, h = T>k B T(f)p g (r)l2yim p , the kinetic energy density profile E kl 



„u~/2 and the turbulent one E u , 



p g u /2; 



the turbulent velocity field u(x) is defined according to the procedures discussed in sect. 6.2. The profiles are distinguished by the corresponding 
labels indicated in the bottom right panel and have been rescaled in units of p^cr;. For the sake of clarity, in each panel only profiles of simulations 
with two different AV parameters (AV and AV 4 ) are shown. 



Vikhlinin & Kravtsov 2007; Piffaretti & Valdarnini 2008; Lau, 
Kravtsov & Nagai 2009). Moreover, the profiles of P lT /(P a -+P,i l ) 
show a weak dependency on the numerical viscosity scheme 
used in the simulations. Differences between the profiles are sig- 
nificant in the cluster cores, but at r — r^o the profiles of runs 
AVo differ by a few per cent from those of runs AV4, indicating 
therefore that pressure terms due to random gas motions are al- 
ready estimated with good approximation in SPH simulations in 
which a standard AV scheme is implemented. The physical con- 
sequences of this result will be discussed in the conclusions, here 
it is worth discussing the origin of the different radial behavior 
of the ratio Pa- 1 {Pa- + Pth) with respect to that of E tU rb/E tot . 
The rms gas velocity cr(r) is defined as 



<r(r) = -y /< « 2 ( r ) > - < u(r) > 2 , 



(34) 



where the field velocity u(r) is calculated according to the SPH 
prescription at the grid points of the spherical shell and brackets 
denote averages performed over the set of points. This definition 
is in accordance with the procedures used to compute the energy 
density profiles; however the gas velocity dispersion calculated 
following the standard definition, i.e. by averaging over the gas 
particles contained within the shell of radial coordinate r and 



thickness Ar, does not differ significantly from that obtained us- 
ing Eq. (34). The different radial behavior of the rms gas velocity 
dispersion <j(r) from that of the local velocity field u(r) can be 
ascribed to the different procedures adopted for computing the 
two velocities, which in turn define the two fields. 

Because of the filtering procedures, the latter is dominated by 
the small-scale random gas motion with a coherence length scale 
~ 200 - 300 kpc. For the rms gas velocity cr(r), the contribution 
of laminar infall motion is not removed and it is mixed with that 
of turbulent motion. Therefore, while the local velocity field u(r) 
can be properly defined as a turbulent velocity field (Vazza et al. 
2009), the dispersion <x(r) is just a measure of the gas velocity 
variance at the scale r. This is the reason why the ratio E tur t,/E tot 
does not show a well-defined dependency with radius, while cr(r) 
increases with r as the streaming motion due to accretion from 
the infalling material becomes progressively more important. A 
similar viewpoint is discussed by Zhu, Feng & Fang (2010, see 
also Maier et al. 2009), who studied the generation of vortic- 
ity fields by means of of hydrodynamic simulations in order to 
measure the level of turbulence in the intergalactic medium. 

To summarize, the notion of turbulent pressure support, 
which accounts for the X-ray mass bias when assuming hydro- 
static equilibrium for the ICM, is somewhat misleading at large 
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cluster radii because the velocity terms which are missed in the 
Jeans equation are dominated by the streaming motion of the 
infalling material, with the local turbulent velocity field progres- 
sively becoming relatively less important in the cluster outskirts. 

6.3. Simulations with cooling and star formation 

The discussion of the previous section was aimed at investigat- 
ing the effects of numerical viscosity on ICM properties when 
the gas evolution of the cluster is driven by a relatively sim- 
ple physics, in which gravity is the sole source of structure for- 
mation. However, a realistic modeling for the physics of the 
ICM must incorporate the possibility for the gas to cool radia- 
tively, turning cold dense clumps of gas into stars. Moreover, 
energy and metal feedback that follows from supernova explo- 
sions should be considered as well. Hydrodynamical simulations 
in which these physical processes have been taken into account 
(Muanwong et al. 2002; Tornatore et al. 2003; Valdarnini 2003; 
Kay et al. 2004; Nagai et al. 2007b) have shown that there is 
substantial agreement between the predicted ICM properties, 
such as the temperature profiles, and the corresponding measure- 
ments, though significant discrepancies still remain (Borgani & 
Kravtsov 2009). In particular, the model overpredicts the star 
mass fraction, which is found in simulations to be a factor ~ 2-3 
higher than that estimated from observations. This is the so- 



called 'overcooling' problem, for which several explanations 
have been proposed (Borgani & Kravtsov 2009), although there 
is not yet a definitive solution. 

An improvement in the AV scheme used in the SPH code is 
expected to significantly modify the ICM properties of the sim- 
ulated cooling clusters. To investigate these effects, the analyses 
of the previous sections are repeated here using the same set 
of cluster samples, but in the simulations the physical modeling 
of the gas now includes radiative cooling and star formation, as 
well as energy and metal feedback from supernovae (Valdarnini 
2006). Our analysis of the results of the simulations will not 
replicate the thorough discussion of Sections 6.1 and 6.2, but 
only the most important results will be presented from the view- 
point of AV when cooling is incorporated into the simulations. 

When applying the spectral analysis of sect. 6.1 to the ve- 
locity fields produced by radiative simulations, a striking feature 
which emerges from the spectral behavior of the transformed 
fields is the presence at small scales of a driving source which 
injects turbulence into the ICM. The density-weighted velocity 
power spectra for the compressive components of the relaxed 
cluster subsample are shown in Fig. 18 and exhibit the property 
of having a spectral behavior which stops decaying at k ~ 10 and 
thereafter is nearly flat at higher wavenumbers, or even grows as 
in the case of cluster 110. These features are in sharp contrast 
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with the corresponding spectra of Fig. 7 for the adiabatic runs 
and are clearly induced by the presence of cooling. 

Previous investigations have considered generation of tur- 
bulence in cluster cores to occur either through the interaction 
of the ambient ICM with buoyantly rising bubbles created from 
AGN jet activities (Churazov et al. 2001; Briiggen & Kaiser 
2002), or induced by orbital motion of galaxies (El-Zant, Kim 
& Kamionkowski 2004; Kim 2007). The origin of turbulence 
is interpreted here as being due to the development of dense 
compact cool gas cores in the inner cluster regions. These are 
characterized by the presence of central gas densities for which 
pipe ~ 10 4 , about a factor ~ 10 higher than in the correspond- 
ing adiabatic runs. Following Fujita, Matsumoto & Wada (2004) 
and Fujita (2005), the interaction of these dense cool cores with 
the bulk gas motion of the surrounding ICM leads then to hy- 
drodynamical instabilities and to the production of turbulence. 
That the origin of the turbulent injection resides in the interac- 
tion of the cluster core with the ambient medium is confirmed 
by the radial dependence of the viscosity parameters (Fig. 19), 
which exhibit a rising ramp in the proximity of the cluster cen- 
ters because of the weak shocks which forms near the core. Due 
to resolution issues, it appears difficult to assess in a quantitative 
way the small-scale behavior of the derived spectra, although it 
seems unlikely that this effect should not be confirmed in simu- 
lations in which the spectra are adequately resolved, as similar 



features are already present in the low resolution runs MR and 
LR. For the same reason no attempt is made here to extract any 
possible dependence of the spectra on the numerical viscosity of 
the simulations, as there is a large scatter at small scales between 
spectra of a given cluster but with different AV settings. 

These results provide strong support however for the plau- 
sibility of models in which radiative cooling in the cluster 
cores is compensated by turbulent heating of the ICM (Fujita, 
Matsumoto & Wada 2004; Fujita 2005). This scenario is obser- 
vationally motivated by the lack of resonant scattering from X- 
ray spectra (Churazov et al. 2004) as well as from the constraints 
on diffusion processes extracted from measured iron abundance 
profiles (Rebusco et al. 2006; David & Nulsen 2008), which 
indicate how in cluster cores the outward diffusion of iron must 
be driven by turbulent gas motion. Heating of the ICM has been 
investigated by Dennis & Chandran (2005), who constructed 
a set of analytical models in which radiative cooling is locally 
balanced by dissipation of turbulent motion, entropy mixing via 
turbulent diffusion and heat conduction. To compensate radiative 
losses, both the heating rates from turbulent diffusion and dissi- 
pation are in general important for establishing thermal equilib- 
rium. The local turbulent heating rate from turbulent dissipation 
is given by 
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Fig. 19. As in Fig. 15, final profiles of the viscosity parameter a are shown for the cooling runs of the perturbed test clusters with different AV 
parameters and resolution. 



where / is the driving length scale and in order to balance ra- 
diative losses estimates derived from simple models (Dennis & 
Chandran 2005; David & Nulsen 2008) require S ~ 100 - 300 
km s -1 . Such a range of values is smaller than that found here 
in the central cluster regions for the local velocity u of the sim- 
ulated clusters by roughly a factor ~ 50%, showing therefore 
that analytical models might overestimate the amount of entropy 
mixing necessary to achieve thermal balance in cluster cores. 
Note, however, that in massive clusters a viable mechanism for 
the dissipation of turbulence is collisionless damping (Brunetti 
& Lazarian 2007), hence for these clusters the contribution of 
the heating rate (35) to the ICM thermal balance might be over- 
estimated. Finally, as indicated by Fig. 19, a certain amount of 
shock heating is present to reduce the effects of cooling in the 
core. 

High levels of turbulence in cluster cores can also have a sig- 
nificant impact on the stability of cold fronts. These structures 
are interpreted as contact discontinuities which originate when 
gas of different entropy is brought into contact (Markevitch & 
Vikhlinin 2007). In the sloshing scenario (ZuHone, Markevitch 
& Johnson 2010, and references therein) this happens because 
of the interaction of the cool core gas with substructures falling 
through the main cluster. Accordingly, cold fronts are created 
as the subsonic sloshing of the central gas causes the encounter 
of gas flows from different directions. However, cold fronts are 



prone to disruption by hydrodynamical instabilities (ZuHone, 
Markevitch & Johnson 2010; Roediger et al. 2010) and it is 
unclear if the level of turbulence found in the cluster cores of 
the simulations shown here can be reconciled with the presence 
of cold fronts as seen in many clusters. Such a issue is beyond 
the scope of this paper and clearly deserves further investiga- 
tions. In particular, hydrodynamic simulations which incorpo- 
rate a proper modeling of the physical viscosity of the ICM can 
be used to investigate the stability of cold fronts in the slosh- 
ing scenario (ZuHone, Markevitch & Johnson 2010). The de- 
rived constraints on the ICM viscosity should then permit to ad- 
dress the viability of turbulent heating models to balance radia- 
tive losses. 

The radial profiles of the energy density ratios, together with 
the relative pressure terms, are shown in Fig. 20. These profiles 
are the analog for the cooling runs of those of Fig. 17 for the 
nonradiative case. A comparison between the two figures shows 
that at r t 0.1r2oo there are no large variations in the E tur b/E tot 
profiles, with the exception of cluster 105. This follows be- 
cause at r ~ 0.1r2oo the steepening in E lur b is accompanied by 
a corresponding steepening in E 10t , so that the ratio E, llr b/E, ot 
retains a profile approximately similar to that of the nonradia- 
tive case. The same behavior is found for the radial profiles of 
Pa-/(Pcr + Pth), which show a very little dependence on the AV 
strength of the simulations at r Z 0.1r2oo- These results hold as 
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Fig. 20. As in Fig. 17, but for the cooling runs. 



well for the relaxed subsample. Because the physical descrip- 
tion of the ICM in these simulations is more realistic than in the 
nonradiative case, these results provide support for the findings 
of the previous section, indicating that numerical viscosity ef- 
fects have little impact on the bias of hydrostatic mass estimates. 
A similar conclusion holds for the entropy profiles, which are 
shown in Fig. 21 and exhibit very little scatter between runs with 
different AV settings. This is interpreted as a consequence of 
galaxy formation processes so that the small amount of entropy 
mixing, which was found in the cluster cores for the adiabatic 
runs, has now been removed. 

The radial profiles of the gas mass fraction f g (< r) = M gas (< 
r)/M c i(< r) are shown in Fig. 22 for the perturbed test clusters. 
The quantity f g (< r) is an increasing function of the radius r, as 
the amount of cooled gas subject to star formation is a decreas- 
ing function of radius because of the dependence of the cooling 
function on the gas density. At r = r^oo ~ 0.6r2oo the value of 
fgas is practically independent of the AV scheme used in the sim- 
ulations, whereas at inner radii it is difficult to extract meaning- 
ful informations about the dependence of f g (< r) on the imple- 
mented AV scheme. For instance, at radii r i 0.1r2oo the highest 
values of / g (< r) are those of runs AV5 for clusters 001, 005 and 
105. The lowest values of f gas are found at r t 0.1r2oo for the 
standard viscosity runs AVq in the case of clusters 005 and 105. 

According the new AV formulation, turbulent heating of the 
ICM is higher because of the suppression of viscous damping of 



random gas motion. This in turn implies a reduced star formation 
efficiency because of a higher level of turbulent pressure, and at 
a given radius the quantity f gas would exhibit higher values for 
those runs where the AV is reduced. The scatter between the pro- 
files of Fig. 22 does not allow to reach firm conclusions about 
this dependency, indicating that a very large simulated sample is 
needed in order to clarify this issue. Finding a statistically sig- 
nificant dependence of f gas on the AV scheme used in the simu- 
lations would support the viewpoint that turbulence might play a 
significant role in solving the overcooling problem, as suggested 
by Zhu, Feng & Fang (2010). 



To summarize, the role of numerical viscosity in simula- 
tions which incorporate radiative cooling is significant in cluster 
cores, where the strength of turbulence driven by hydrodynam- 
ical instabilities is strongly amplified with respect to the non- 
radiative case, but is negligible in establishing the shape of the 
ICM profiles at cluster radii r i 0.1 ^oo- Results from the pre- 
vious section support the view that simulations with a much 
larger number of gas particles ( > 256 3 ) are necessary in order 
to consistently investigate the first of these issues, while at radii 
r i 0.1r2oo the stability of the ICM profiles presented here ap- 
pears quite robust. 
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7. Summary and Conclusions 

In this paper results from hydrodynamical SPH simulations have 
been presented aimed at investigating the role of the SPH nu- 
merical viscosity scheme in determining ICM properties of sim- 
ulated galaxy clusters. For doing this, the standard AV formu- 
lation of SPH has been modified according to the proposal of 
Morris & Monaghan (1997), where each particle has its own 
viscosity parameter a{t) which evolves in time under the local 
shock conditions. This time-dependent AV scheme has the main 
advantage of properly reducing the amount of AV present in the 
fluid in regions far from shocks, with respect the standard AV 
scheme where the corresponding parameters are held constant 
throughout the simulation domain. The new time-dependent AV 
scheme has been implemented in an SPH code and several nu- 
merical tests have been performed to assess the validity of the 
code and its shock resolution capabilities. The chosen tests were 
the Riemann shock tube problem and 3D collapse of a cold gas 
sphere. Both of these problems have been widely studied in the 
literature and results from previous simulations provide refer- 
ence profiles against which to compare the profiles extracted 
from the simulations produced here. 

For the tests considered the results of simulations performed 
using the standard AV scheme are in good agreement with previ- 
ous findings, showing the validity of the code. In order to inves- 
tigate the impact of the time-dependent AV scheme in reducing 
numerical viscosity effects, for a given test case a set of SPH 



simulations which incorporate the new AV formulation has been 
constructed, in which the runs use the same initial conditions 
but differ in the settings of the AV parameters which control the 
time evolution of the individual viscosity coefficients. The re- 
sults of these simulations show profiles which are very similar 
to those produced by the fixed AV runs, indicating that the new 
AV scheme has shock resolution capabilities similar to those of 
the standard one and at the same time a reduced AV far from 
the shocks. Moreover, the time evolution of the viscosity param- 
eter a(t) is consistent with theoretical expectations and at the 
shock front the coefficients reach peak values which are nearly 
independent of the chosen AV settings. Finally, in the cold gas 
sphere, the entropy profiles of the time-dependent AV runs ex- 
hibit reduced pre-shock heating at the shock front and better 
agreement with the reference ID PPM numerical solution. 

Having established the reliability of the code and the effec- 
tiveness of the new AV formulation, the impact of AV on ICM 
properties of simulated clusters has been investigated by extract- 
ing results from a large ensemble of hydrodynamical cluster 
simulations. The ensemble has been constructed by collecting 
samples of simulated clusters, each of the samples being real- 
ized by performing simulations for the same set of cluster ini- 
tial conditions but using a different setting of the AV parame- 
ters in the SPH code. The baseline cluster sample with which 
the simulation ensemble has been constructed consisted of eight 
clusters chosen from a cosmological simulation ensemble with 
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cluster virial masses ranging from M200 ~ 5 • 1O 13 /i~'M up to 
M200 ~ 6 • 10 14 /z -1 M G , the selection criterion being that of con- 
structing a representative sample of cluster dynamical states and 
masses. The stability of the results against the numerical resolu- 
tion of the simulations has been tested by constructing a set of 
mirror ensembles in which the only difference in the cluster sim- 
ulations with respect to the corresponding one of the reference 
ensemble was in the number of simulation particles. Moreover, 
the ensembles of adiabatic cluster simulations have been repli- 
cated following the same prescriptions but with the gas being 
allowed to cool radiatively and being subject to star formation. 
The final outcome of these procedures is a set of hydrodynamical 
cluster simulations which allow the dependence of ICM gas ve- 
locities on the AV scheme used in the simulations to be studied 
in a systematic way. Here, the main results are summarized. 

The velocity Fourier spectra of the simulated clusters have 
in common certain features which allow several conclusions 
to be drawn about the generation of random gas velocities. 
All of the spectra exhibit a maximum at the smallest spec- 
tral cube wavenumber, showing that merging and accretion pro- 
cesses driven by gravity at cluster scales are the primary sources 
of energy injection into the ICM. The spectra begin to man- 
ifest their dependence on dissipative effects at length scales 
Idiss ~ ^200/10 ~ 100 - 300 kpc, below which the less dissi- 
pative spectra are those extracted from runs in which the AV 



settings correspond to the shortest decay time scale for the vis- 
cosity parameter a( f). The range of values for l^ ss appears to be 
in accordance with previous findings (Subramanian, Shukurov 
& Haugen 2006; Maier et al. 2009), suggesting that the de- 
tected spectral features are not an artifact due to resolution ef- 
fects. These behaviors are consistent with those derived from 
the structure function radial profiles and are common to all of 
the clusters, regardless of the dynamical state of the cluster. 

Spectral decomposition into a transverse and a longitudinal 
part allows the resolution dependence of the spectra to be stud- 
ied and useful insights to be obtained into the resolution require- 
ments which must be fulfilled in order to adequately resolve tur- 
bulence in SPH simulations of galaxy clusters. The longitudinal 
component of the velocity power spectrum E c (k) is found to ex- 
hibit an approximate Kolgomorov scaling over nearly a decade, 
whilst the shearing part E s (k) has a much steeper fall off. In anal- 
ogy with previous findings (Federrath et al. 2010), this behavior 
is interpreted as a resolution effect, as the number of computa- 
tional elements needed to adequately resolve rotational motion 
is higher than that necessary for longitudinal modes, where only 
one degree of freedom is present. Using scaling arguments de- 
rived from previous simulations (K09, Price & Federrath 2010) 
a number of N ~ 512 3 gas particles is then needed to adequately 
describe turbulence driven spectra in SPH simulated clusters. 
However, this estimate might be too pessimistic, as the SPH res- 
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olution is not constant but increases with the particle density to- 
ward the cluster center, although it appears unlikely that less than 
~ 256 3 particles can be used to describe E s ik) over a decade. 

Entropy is a fundamental thermodynamic variable and the 
corresponding average radial profiles have been chosen to in- 
vestigate the dependence of the final ICM thermal properties on 
the AV scheme used in the simulations. As a function of the AV 
simulation parameters for the entropy profiles of a given clus- 
ter, a very small scatter is found with the largest differences in 
most cases being £ 50% within radial distances riO.l r2oo ■ The 
weak dependency of the entropy profiles on the AV parameters 
of the simulations indicates that the amount of entropy mixing in 
cluster cores due to numerical viscosity effects can be considered 
negligible. As a consequence, the discrepancy between SPH and 
Eulerian codes in the core entropy produced in adiabatic clus- 
ter simulations can be primarily related to the treatment of fluid 
discontinuities in the SPH formulation, rather than to AV effects. 

At variance with the entropy profiles, the radial behavior of 
the turbulent energy density profile E, ur b(r) exhibits a signifi- 
cant dependence on the chosen AV scheme. There is an increase 
in the turbulent-to-total energy density ratio as less dissipative 
AV runs are considered, with the largest variations occurring 
at r < 0.1r2oo- The range of values for the ratio E tur \,jE tot is 
in accordance with previous findings obtained using the AMR 
Eulerian code ENZO (Vazza et al. 2009) and provides sup- 
port for the view that differences between results extracted from 
SPH and grid based codes can be reconciled using an adequate 
number of simulation particles in SPH simulations and properly 
treating the AV. 

A striking result is the lack of any significant dependence at 
large cluster radii of the rms gas velocity cr(r) on the AV param- 
eters used in the simulations. As this quantity dominates in the 
Jeans equation the mass correction terms to the hydrostatic equi- 
librium equation, it follows that numerical viscosity effects can 
be considered negligible in determining the corrections to the 
hydrostatic X-ray cluster mass estimates. Hence, previous SPH 
simulations in which the AV is treated according to the standard 
formulation, already provide numerically reliable estimates for 
the X-ray mass bias (Piffaretti & Valdarnini 2008). One impor- 
tant consequence of this result is that the agreement between 
the measured hydrostatic X-ray and weak lensing mass ratio 
(Mahdavi et al. 2008; Zhang et al. 2008, 2010) with the X-ray 
mass bias predicted by simulations (Zhang et al. 2010), strongly 
supports the already outlined scenario (Piffaretti & Valdarnini 
2008), in which the level of non-thermal pressure present in the 
ICM is not significant and the gas physics outside cluster cores 
is well described by the current SPH simulations. 

In simulations which incorporate radiative cooling, a signifi- 
cant feature is the presence in the inner cluster regions of a much 
higher level of turbulence, which is produced by the hydrody- 
namical instabilities generated by the interaction of the compact 
cool core with the ambient medium. Moreover, the range of val- 
ues of the local turbulent velocity is in accordance with the con- 
straints required by the turbulent heating model in order to bal- 
ance radiative losses (Dennis & Chandran 2005). 

To address in a quantitative way the viability of the tur- 
bulent heating model, a proper treatment would be needed of 
the small-scale fluid instabilities which develop near the cluster 
cores. This requires the use of simulations with a much larger 
number of particles and at the same time incorporating into the 
SPH equations the thermal diffusion terms which have been pro- 
posed (Price 2008; Wadsley, Veeravalli & Couchman 2008) to 
correctly model fluid instabilities in SPH. The latter is particu- 
larly relevant in order to reliably estimate the radial profile of the 



local velocity S(r), which is an important factor in establishing 
the local balance between the turbulent heating rate and radia- 
tive losses (Dennis & Chandran 2005). Such simulations will 
still miss some basic ingredients of the physics describing the 
thermodynamics of the ICM in cluster cores, in particular the 
effects due to thermal conduction, heating from AGN jets, phys- 
ical viscosity and instabilities driven by the presence of magnetic 
fields (Parrish, Quataert & Sharma 2010). Nonetheless, the in- 
formation provided by these simulations will probably shed light 
on the role played by turbulence in establishing the observed bi- 
modal distribution for cluster core entropies (Pratt et al. 2010) 
and the morphological distinction between cool-core and non 
cool-core clusters. 

To summarize, the results presented here indicate that turbu- 
lence is likely to play a significant role in some open issues of 
ICM physics, such as the cooling flow and the overcooling prob- 
lems. To properly investigate these issues it would be necessary 
however to further improve the hydrodynamic SPH equations, 
so as to correctly treat fluid discontinuities, and to use a much 
larger number of particles. For a given test cluster, a comparison 
between the results produced by these simulations and those ex- 
tracted from simulations obtained using an Eulerian AMR code 
is likely to give interesting insights into both the physical prob- 
lem and also the numerical framework. 
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